Setting Directory

We first set the directory for the EPA plots. Installing required packages.

Installing Packages

Loading Libraries

We load all required libraries together.

library(sf)
library(corrr)
## Warning: package 'corrr' was built under R version 4.1.1
library(tmap)
library(rgdal)
library(dplyr)
library(raster)
library(data.table)
library(sp)
library(tidyverse)
library(tidycensus)
library(psych)
library(tigris)
library(leaflet)
library(maps)
library(gstat)
library(chron)
library(maptools)
library(ggcorrplot)
library(ggpubr)
library(Hmisc)
library(corrplot)

##Inputting CSVs

We inputting the various csv files of different pollutants. (What are they, where they came from, what all they did, what’s the goal, etc.). (Show fancier tables, like cable)

#read in and wrangle epa data for SO2 
so2_epa <- read.csv("SO2.csv")
so2_epa_df <- so2_epa %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(SO2,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("SO2")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


so2_epa_df <- so2_epa_df %>% dplyr::filter(COUNTY == "Philadelphia")

#read in and wrangle epa data for PM2.5 
pm_epa <- read.csv("PM25.csv")

pm_epa_df <- pm_epa %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Mean.PM2.5.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("PM2.5")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_epa_df <- pm_epa_df %>% dplyr::filter(COUNTY == "Philadelphia")

#read in and wrangle epa data for Ozone 
pm_oz <- read.csv("Ozone.csv")

pm_oz_df <- pm_oz %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Max.8.hour.Ozone.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("OZ")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_oz_df <- pm_oz_df %>% dplyr::filter(COUNTY == "Philadelphia")


#read in and wrangle epa data for NO2 
pm_no2 <- read.csv("NO2.csv")

pm_no2_df <- pm_no2 %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Max.1.hour.NO2.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("NO2")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_no2_df <- pm_no2_df %>% dplyr::filter(COUNTY == "Philadelphia")


#read in and wrangle epa data for Lead 
pm_pb <- read.csv("Pb.csv")

pm_pb_df <- pm_pb %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Mean.Pb.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("Pb")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_pb_df <- pm_pb_df %>% dplyr::filter(COUNTY == "Philadelphia")

#read in and wrangle epa data for CO 
pm_co <- read.csv("CO.csv")

pm_co_df <- pm_co %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Max.8.hour.CO.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("CO")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_co_df <- pm_co_df %>% dplyr::filter(COUNTY == "Philadelphia")


#read in and wrangle epa data for PM10 
pm_10 <- read.csv("PM10.csv")

pm_10_df <- pm_10 %>% 
  dplyr::mutate(
    Date = as.factor(gsub("/00","/",strftime(as.Date(Date, format="%m/%d/%Y"), format="%m/%d/%Y"))),
    Day = gsub("/([^/]*)$","",Date),
    Year = as.factor(paste0("20",gsub(".*/","",Date))),
    Mean_value = round(Daily.Mean.PM10.Concentration,2),
    SiteName = as.factor(gsub(".*[(]","(","Site Name")),
    Parameter_name = as.character("PM10")) %>% dplyr::rename(Latitude = SITE_LATITUDE, Longitude = SITE_LONGITUDE)


pm_10_df <- pm_10_df %>% dplyr::filter(COUNTY == "Philadelphia")

head(pm_10_df)
##         Date Source   Site.ID POC Daily.Mean.PM10.Concentration    UNITS
## 1 12/03/2020    AQS 421010048   2                            15 ug/m3 SC
## 2 12/04/2020    AQS 421010048   2                            21 ug/m3 SC
## 3 12/05/2020    AQS 421010048   2                            13 ug/m3 SC
## 4 12/06/2020    AQS 421010048   2                             8 ug/m3 SC
## 5 12/07/2020    AQS 421010048   2                            12 ug/m3 SC
## 6 12/08/2020    AQS 421010048   2                            15 ug/m3 SC
##   DAILY_AQI_VALUE              Site.Name DAILY_OBS_COUNT PERCENT_COMPLETE
## 1              14 North East Waste (NEW)               1              100
## 2              19 North East Waste (NEW)               1              100
## 3              12 North East Waste (NEW)               1              100
## 4               7 North East Waste (NEW)               1              100
## 5              11 North East Waste (NEW)               1              100
## 6              14 North East Waste (NEW)               1              100
##   AQS_PARAMETER_CODE    AQS_PARAMETER_DESC CBSA_CODE
## 1              81102 PM10 Total 0-10um STP     37980
## 2              81102 PM10 Total 0-10um STP     37980
## 3              81102 PM10 Total 0-10um STP     37980
## 4              81102 PM10 Total 0-10um STP     37980
## 5              81102 PM10 Total 0-10um STP     37980
## 6              81102 PM10 Total 0-10um STP     37980
##                                     CBSA_NAME STATE_CODE        STATE
## 1 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
## 2 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
## 3 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
## 4 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
## 5 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
## 6 Philadelphia-Camden-Wilmington, PA-NJ-DE-MD         42 Pennsylvania
##   COUNTY_CODE       COUNTY Latitude Longitude   Day   Year Mean_value  SiteName
## 1         101 Philadelphia 39.99139 -75.08083 12/03 202020         15 Site Name
## 2         101 Philadelphia 39.99139 -75.08083 12/04 202020         21 Site Name
## 3         101 Philadelphia 39.99139 -75.08083 12/05 202020         13 Site Name
## 4         101 Philadelphia 39.99139 -75.08083 12/06 202020          8 Site Name
## 5         101 Philadelphia 39.99139 -75.08083 12/07 202020         12 Site Name
## 6         101 Philadelphia 39.99139 -75.08083 12/08 202020         15 Site Name
##   Parameter_name
## 1           PM10
## 2           PM10
## 3           PM10
## 4           PM10
## 5           PM10
## 6           PM10

Creating Base Raster

## 0.1km Rasters for Philadelphia only
getRaster_ph <- function(data, par_name){
  ## filter data to contiguous US
  ## get one row per location by averaging all entries
  dat <- dplyr::filter(data,
                       Parameter_name == par_name) %>%
    dplyr::filter(!is.na(Mean_value), Mean_value >= 0) %>%
    dplyr::group_by(Latitude, Longitude) %>%
    dplyr::summarise(avg = mean(Mean_value, na.rm = TRUE))
  
  if(nrow(dat) > 0){
    coordinates(dat) = ~Longitude+Latitude
    crs(dat) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
    
    ## make base raster
    #find distance between the latitude and longitudes and convert to km (*111 for 1km)
    r <- raster(nrow = 301, ncol = 360, extent(-75.28027, -74.95576, 39.867, 40.13799)) #0.1km
    crs(r) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
    ## generate raster (idw with 5 nearest sites)
    gs <- gstat(formula=avg~1, data=dat, nmax = 5)
    nn <- interpolate(r, gs)
    
    return(nn)
  } else{
    print("zero")
    #find distance between the latitude and longitudes and convert to km (*111 for 1km)
    r <- raster(nrow = 301, ncol = 360, extent(-75.28027, -74.95576, 39.867, 40.13799))
    crs(r) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
    values(r) <- rep(-1, 108360)
    return(r)
  }
}

Writing Rasters

#PH shape file
ph_shapefile <- readOGR("./Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/karansampath/Desktop/Research/PURM/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.shp", layer: "Neighborhoods_Philadelphia"
## with 158 features
## It has 5 fields
ph <- spTransform(x = ph_shapefile, CRSobj = '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')

#Census Shapefile
census_shapefile <- readOGR("./census_shapefile/tl_2020_42_tract.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/karansampath/Desktop/Research/PURM/census_shapefile/tl_2020_42_tract.shp", layer: "tl_2020_42_tract"
## with 3446 features
## It has 12 fields
## Integer64 fields read as strings:  ALAND AWATER
#make pm.25 raster
pm_ph_daily_brick <- brick(getRaster_ph(pm_epa_df,"PM2.5"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
pm_ph_daily <- getRaster_ph(pm_epa_df,"PM2.5")
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
pm_ph_daily_brick[pm_ph_daily_brick == -1] <- NA
names(pm_ph_daily_brick) <- c("pm2.5")
c <- crop(pm_ph_daily_brick, extent(ph))
pm25_raster <- raster::mask(c, ph) 
writeRaster(pm25_raster,"ph_pm_daily_0.1km_brick.tif",overwrite=TRUE)

#make so2 raster
SO2_ph_daily_brick <- brick(getRaster_ph(so2_epa_df,"SO2"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
SO2_ph_daily_brick[SO2_ph_daily_brick == -1] <- NA
names(SO2_ph_daily_brick) <- c("so2")
#crop by shape
c <- crop(SO2_ph_daily_brick, extent(ph))
so2_raster <- raster::mask(c, ph) 
writeRaster(so2_raster,"ph_SO2_daily_0.1km_brick.tif",overwrite=TRUE)


#make CO raster
CO_ph_daily_brick <- brick(getRaster_ph(pm_co_df,"CO"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
CO_ph_daily_brick[CO_ph_daily_brick == -1] <- NA
names(CO_ph_daily_brick) <- c("CO")
#crop by shape
c <- crop(CO_ph_daily_brick, extent(ph))
co_raster <- raster::mask(c, ph) 
writeRaster(co_raster,"ph_CO_daily_0.1km_brick.tif",overwrite=TRUE)

#make Ozone raster
oz_ph_daily_brick <- brick(getRaster_ph(pm_oz_df,"OZ"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
oz_ph_daily_brick[oz_ph_daily_brick == -1] <- NA
names(oz_ph_daily_brick) <- c("OZ")
#crop by shape
c <- crop(oz_ph_daily_brick, extent(ph))
o3_raster <- raster::mask(c, ph) 
writeRaster(o3_raster,"ph_oz_daily_0.1km_brick.tif",overwrite=TRUE)

#make NO2 raster
NO2_ph_daily_brick <- brick(getRaster_ph(pm_no2_df,"NO2"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
NO2_ph_daily_brick[NO2_ph_daily_brick == -1] <- NA
names(NO2_ph_daily_brick) <- c("NO2")
#crop by shape
c <- crop(NO2_ph_daily_brick, extent(ph))
no2_raster <- raster::mask(c, ph) 
writeRaster(no2_raster,"ph_NO2_daily_0.1km_brick.tif",overwrite=TRUE)



#make Pb raster
pb_ph_daily_brick <- brick(getRaster_ph(pm_pb_df,"Pb"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
pb_ph_daily_brick[pb_ph_daily_brick == -1] <- NA
names(pb_ph_daily_brick) <- c("Pb")
#crop by shape
c <- crop(pb_ph_daily_brick, extent(ph))
pb_raster <- raster::mask(c, ph) 
writeRaster(pb_raster,"ph_pb_daily_0.1km_brick.tif",overwrite=TRUE)


#make PM10 raster
pm10_ph_daily_brick <- brick(getRaster_ph(pm_10_df,"PM10"))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
## [inverse distance weighted interpolation]
pm10_ph_daily_brick[pm10_ph_daily_brick == -1] <- NA
names(pm10_ph_daily_brick) <- c("PM10")
#crop by shape
c <- crop(pm10_ph_daily_brick, extent(ph))
pm10_raster <- raster::mask(c, ph) 
writeRaster(pm10_raster,"ph_pm10_daily_0.1km_brick.tif",overwrite=TRUE)
pm25_poly <- raster::extract(pm25_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
pm10_poly <- raster::extract(pm10_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
no2_poly <- raster::extract(no2_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
o3_poly <- raster::extract(o3_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
co_poly <- raster::extract(co_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
so2_poly <- raster::extract(so2_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
pb_poly <- raster::extract(pb_raster, census_shapefile)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
# Note output is a list of 286 vectors, corresponding to the 286 polygons in the sepa zip/muni shapefile. Each vector contains all the raster values falling within a polygon.
str(census_shapefile, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  3446 obs. of  12 variables:
##   ..@ polygons   :List of 3446
##   ..@ plotOrder  : int [1:3446] 1109 1557 42 13 301 482 3105 1544 481 372 ...
##   ..@ bbox       : num [1:2, 1:2] -80.5 39.7 -74.7 42.5
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   ..$ comment: chr "TRUE"
str(pm25_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 8.09 8.09 8.09 8.09 8.09 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 8.07 8.08 8.08 8.08 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:31, 1] 7.91 7.9 7.9 7.89 7.88 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:134, 1] 7.85 7.85 7.86 7.85 7.85 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:20, 1] 7.62 7.62 7.62 7.62 7.62 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:46, 1] 7.62 7.61 7.61 7.61 7.61 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:47, 1] 7.57 7.58 7.58 7.58 7.59 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 7.86 7.86 7.86 7.86 7.86 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:199, 1] 7.87 7.87 7.87 7.87 7.87 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:37, 1] 7.49 7.49 7.49 7.49 7.49 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:107, 1] 7.69 7.7 7.7 7.68 7.68 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:22, 1] 7.61 7.61 7.61 7.61 7.61 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:65, 1] 7.84 7.84 7.84 7.84 7.84 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : num [1:225, 1] 7.85 7.85 7.85 7.85 7.85 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "pm2.5"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(pm10_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:31, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:134, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:20, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:46, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:47, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:199, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:37, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:107, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:22, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:65, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : num [1:225, 1] 19 19 19 19 19 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "PM10"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(no2_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 22.6 22.6 22.6 22.6 22.6 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 22.6 22.6 22.6 22.6 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:31, 1] 22.5 22.5 22.5 22.5 22.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:134, 1] 23.5 23.5 23.6 23.5 23.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:20, 1] 22.4 22.4 22.4 22.4 22.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:46, 1] 22.5 22.5 22.5 22.5 22.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:47, 1] 22.4 22.4 22.4 22.4 22.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 22.4 22.4 22.4 22.4 22.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:199, 1] 24 24 24 24 24 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:37, 1] 22.5 22.5 22.5 22.5 22.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:107, 1] 22.1 22.1 22.1 22.1 22.1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:22, 1] 22.4 22.4 22.4 22.4 22.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:65, 1] 22.4 22.4 22.4 22.4 22.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : num [1:225, 1] 23.6 23.6 23.6 23.6 23.6 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "NO2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(o3_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 0.0342 0.0342 0.0342 0.0342 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 0.0342 0.0342 0.0342 0.0342 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:31, 1] 0.0342 0.0342 0.0342 0.0342 0.0342 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:134, 1] 0.0363 0.0363 0.0363 0.0362 0.0363 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:20, 1] 0.0343 0.0343 0.0343 0.0343 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:46, 1] 0.0343 0.0343 0.0343 0.0343 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:47, 1] 0.0343 0.0343 0.0343 0.0343 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 0.034 0.034 0.034 0.034 0.034 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:199, 1] 0.0365 0.0365 0.0366 0.0365 0.0365 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:37, 1] 0.0343 0.0343 0.0343 0.0343 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:107, 1] 0.0342 0.0342 0.0342 0.0342 0.0342 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:22, 1] 0.0343 0.0343 0.0343 0.0343 0.0343 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:65, 1] 0.034 0.034 0.034 0.034 0.034 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : num [1:225, 1] 0.036 0.036 0.036 0.036 0.036 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "OZ"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(co_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 0.603 0.603 0.603 0.603 0.604 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 0.6 0.601 0.601 0.601 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:31, 1] 0.592 0.591 0.59 0.588 0.587 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:134, 1] 0.539 0.54 0.54 0.538 0.539 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:20, 1] 0.523 0.522 0.525 0.523 0.521 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:46, 1] 0.563 0.563 0.562 0.56 0.559 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:47, 1] 0.504 0.503 0.501 0.5 0.499 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 0.503 0.501 0.499 0.498 0.496 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:199, 1] 0.574 0.574 0.574 0.574 0.575 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:37, 1] 0.556 0.555 0.554 0.554 0.553 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:107, 1] 0.421 0.419 0.416 0.429 0.427 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:22, 1] 0.532 0.531 0.529 0.532 0.53 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:65, 1] 0.472 0.471 0.47 0.469 0.469 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : num [1:225, 1] 0.547 0.547 0.545 0.547 0.547 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "CO"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(so2_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 1.26 1.27 1.27 1.26 1.27 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 1.26 1.26 1.26 1.26 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:31, 1] 1.3 1.3 1.3 1.3 1.3 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:134, 1] 1.16 1.16 1.16 1.16 1.16 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:20, 1] 1.33 1.33 1.33 1.33 1.33 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:46, 1] 1.35 1.35 1.35 1.35 1.35 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:47, 1] 1.34 1.34 1.34 1.34 1.34 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 1.18 1.18 1.17 1.17 1.17 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:199, 1] 1.16 1.16 1.16 1.16 1.16 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:37, 1] 1.37 1.37 1.37 1.37 1.36 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:107, 1] 1.16 1.16 1.16 1.17 1.17 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:22, 1] 1.34 1.34 1.34 1.34 1.34 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:65, 1] 1.16 1.16 1.16 1.16 1.16 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : num [1:225, 1] 1.18 1.18 1.18 1.18 1.18 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "so2"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
str(pb_poly)
## List of 3446
##  $ : NULL
##  $ : num [1:121, 1] 0.00272 0.00271 0.0027 0.0027 0.0027 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:28, 1] NA 0.00277 0.00275 0.00275 0.00275 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:31, 1] 0.00362 0.00366 0.00371 0.00375 0.00379 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:134, 1] 0.00507 0.00508 0.00508 0.00507 0.00507 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:20, 1] 0.00732 0.0073 0.00742 0.0074 0.00739 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:46, 1] 0.0073 0.00743 0.00744 0.00745 0.00746 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:47, 1] 0.0082 0.00817 0.00814 0.00812 0.00809 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:532, 1] NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : num [1:64, 1] 0.00398 0.00399 0.00401 0.00402 0.00404 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:199, 1] 0.0053 0.0053 0.00529 0.0053 0.0053 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:37, 1] 0.0085 0.00851 0.00851 0.00852 0.00852 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:107, 1] 0.00562 0.00563 0.00564 0.00563 0.00564 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:22, 1] 0.00757 0.00756 0.00755 0.00766 0.00765 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:65, 1] 0.00437 0.00438 0.00439 0.0044 0.00442 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : num [1:225, 1] 0.00523 0.00523 0.0052 0.00523 0.00523 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "Pb"
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##   [list output truncated]
# Find the mean for each polygon
pm25_means <- unlist(lapply(pm25_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
pm10_means <- unlist(lapply(pm10_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
no2_means <- unlist(lapply(no2_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
o3_means <- unlist(lapply(o3_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
co_means <- unlist(lapply(co_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
so2_means <- unlist(lapply(so2_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
pb_means <- unlist(lapply(pb_poly, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))


# Append mean pm to SEPA data
census_shapefile@data$pm25 <- pm25_means
census_shapefile@data$pm10 <- pm10_means
census_shapefile@data$no2 <- no2_means
census_shapefile@data$o3 <- o3_means
census_shapefile@data$co <- co_means
census_shapefile@data$so2 <- so2_means
census_shapefile@data$pb <- pb_means
df_census <- census_shapefile@data
census_final <- df_census %>% dplyr::select(GEOID, pm25, pm10, no2, o3, co, so2, pb) %>% drop_na()
head(census_final)
##         GEOID     pm25     pm10      no2         o3        co      so2
## 1 42101012204 8.117772 18.96552 22.57388 0.03423595 0.6066132 1.271117
## 2 42101012203 8.085777 18.96552 22.56867 0.03422197 0.6020866 1.259511
## 3 42101013602 7.849828 18.96552 22.52960 0.03419239 0.5835414 1.298763
## 4 42101034502 7.857126 18.96552 23.60571 0.03635065 0.5430756 1.158944
## 5 42101000902 7.621876 18.96552 22.39748 0.03428038 0.5180747 1.326588
## 6 42101001201 7.592708 18.96552 22.48680 0.03428508 0.5561374 1.353303
##            pb
## 1 0.002596183
## 2 0.002719788
## 3 0.003989176
## 4 0.005100759
## 5 0.007401185
## 6 0.007640434
myLabelFormat = function(prefix = "", suffix = "", between = " &ndash; ", digits = 3, 
                         big.mark = ",", transform = identity, t.val = Inf) {
  formatNum <- function(x) {
    format(round(transform(x), digits), trim = TRUE, scientific = FALSE, 
           big.mark = big.mark)
  }
  function(type, ...) {
    switch(type, numeric = (function(cuts) {
      cuts <- sort(cuts, decreasing = T) #just added
      paste0(prefix, formatNum(cuts), ifelse(cuts == t.val, "+", ""))
    })(...), bin = (function(cuts) {
      n <- length(cuts)
      paste0(prefix, formatNum(cuts[-n]), between, formatNum(cuts[-1]), 
             suffix)
    })(...), quantile = (function(cuts, p) {
      n <- length(cuts)
      p <- paste0(round(p * 100), "%")
      cuts <- paste0(formatNum(cuts[-n]), between, formatNum(cuts[-1]))
      paste0("<span title=\"", cuts, "\">", prefix, p[-n], 
             between, p[-1], suffix, "</span>")
    })(...), factor = (function(cuts) {
      paste0(prefix, as.character(transform(cuts)), suffix)
    })(...))
  }
}

Plotting Maps

make_map <- function(ras.t,title="PM2.5",limt,all_values,monitors_ph,palette){
  if (palette %in% c("heat","HEAT","h","H")){
    palette1 <- colorNumeric(rev(heat.colors(8)), all_values, na.color = "transparent") #terrain
    palette2 <- colorNumeric(heat.colors(8), all_values, na.color = "transparent")
  } else if (palette %in% c("terrain","TERRAIN","t","T")) {
    palette1 <- colorNumeric(rev(terrain.colors(8)), all_values, na.color = "transparent")
    palette2 <- colorNumeric(terrain.colors(8), all_values, na.color = "transparent")
  }
map1<-leaflet(ph) %>% 
    setView(-75.14, 40, zoom = 11) %>%
    addRasterImage(x = ras.t, colors = palette1, method = "ngb", opacity = 0.7) %>%
    addPolygons(color = "black", weight = 1, fillColor = "transparent") %>%
    addCircleMarkers(lng = monitors_ph$Longitude, lat = monitors_ph$Latitude, radius = 5,label = monitors_ph$SiteName) %>%
    addLegend(pal = palette2,
              values = all_values,
              labFormat = myLabelFormat(t.val = limt), title = title,
              position = "bottomleft")
map1
}

Making the Bricks

pm2.5.brick <- brick("ph_pm_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pm2.5.brick[[1]], c(11, Inf, 11)) #8
pm25_data <- values(pm.t1)
epa_sites <- pm_epa_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
epa_main <- pm_epa_df %>% dplyr::group_by(Latitude, Longitude) %>% dplyr::summarise(m = mean(Mean_value))
## `summarise()` has grouped output by 'Latitude'. You can override using the `.groups` argument.
P1 <- make_map(pm.t1,title="PM2.5",11,pm25_data,epa_sites,palette="h") 
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P1
pm2.5.brick <- brick("ph_CO_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t2 <- reclassify(pm2.5.brick[[1]], c(11, Inf, 11)) #8
pm25_data <- values(pm.t1)
epa_sites <- pm_co_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P2 <- make_map(pm.t1,title="CO",11,pm25_data,epa_sites,palette="h") 
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P2
pmno2.brick <- brick("ph_NO2_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pmno2.brick[[1]], c(11, Inf, 11)) #8
pmno2_data <- values(pm.t1)
epa_sites <- pm_no2_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P3 <- make_map(pm.t1,title="NO2",11,pmno2_data,epa_sites,palette="h")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P3
pm2.5.brick <- brick("ph_SO2_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pm2.5.brick[[1]], c(11, Inf, 11)) #8
pm25_data <- values(pm.t1)
epa_sites <- so2_epa_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P3 <- make_map(pm.t1,title="SO2",11,pm25_data,epa_sites,palette="h") #save image as precovid 2020
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P3
pm2.5.brick <- brick("ph_oz_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pm2.5.brick[[1]], c(11, Inf, 11)) #8
pm25_data <- values(pm.t1)
epa_sites <- pm_oz_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P4 <- make_map(pm.t1,title="Ozone",11,pm25_data,epa_sites,palette="h") #save image as precovid 2020
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P4
pm2.5.brick <- brick("ph_pm10_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pm2.5.brick[[1]], c(11, Inf, 11)) #8
pm25_data <- values(pm.t1)
epa_sites <- pm_10_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P5 <- make_map(pm.t1,title="PM10",11,pm25_data,epa_sites,palette="h") #save image as precovid 2020
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P5
pb.brick <- brick("ph_pb_daily_0.1km_brick.tif")
#Use max value to reclassify
pm.t1 <- reclassify(pb.brick[[1]], c(11, Inf, 11)) #8
pb_data <- values(pm.t1)
epa_sites <- pm_pb_df %>% dplyr::select(SiteName, Latitude, Longitude) %>% unique()
P6 <- make_map(pm.t1,title="Pb",11,pb_data,epa_sites,palette="h")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
P6
## Make a table with tracts, automatically should come with GEOIDs once you convert bricks to table (need to find GEOIDs and values for the entire raster)
library(htmltools)
leaflet_grid <- 
  tagList(
    tags$table(width = "100%",
      tags$tr(
          tags$td(P1),
          tags$td(P2)
      ),
      tags$tr(
        tags$td(P3),
        tags$td(P4)
      )
    )
  )
browsable(leaflet_grid)
library(gstat)                                         
library(sp)
z <- stack(pm.t1, pm.t2)
r <- calc(z, fun=function(x) cor(x[1], x[2], method='pearson'))
plot(pm.t1)

Plotting Census Data

#census_api_key("3304af6f650f9f6e42986fbc8497d3b58ac1250e", install=TRUE, overwrite = TRUE)
v17 <- load_variables(dataset = "acs5", year = 2019)
view(v17)
race <- get_acs(geography = "tract", year = 2019, variables = c("B03002_001", "B03002_001","B03002_002","B03002_003","B03002_004","B03002_005", "B03002_006", "B03002_007","B03002_008", "B03002_009", "B03002_010", "B03002_011"), state = "PA", county = c("Philadelphia", "Chester", "Montgomery", "Bucks", "Delaware"))
## Getting data from the 2015-2019 5-year ACS
rtracts <- tracts(state = 'PA', county = c("Philadelphia", "Chester" , "Montgomery", "Bucks", "Delaware"), cb = TRUE)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |======================================================================| 100%
class(race)
## [1] "tbl_df"     "tbl"        "data.frame"
class(rtracts)
## [1] "sf"         "data.frame"
plot(st_geometry(rtracts))

race_merged<- geo_join(rtracts, race, "GEOID", "GEOID")
## Warning: `group_by_()` was deprecated in dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
race_merged
## Simple feature collection with 998 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.13645 ymin: 39.72156 xmax: -74.7215 ymax: 40.60858
## Geodetic CRS:  NAD83
## First 10 features:
##    STATEFP COUNTYFP TRACTCE             AFFGEOID       GEOID  NAME.x LSAD
## 1       42      101  014500 1400000US42101014500 42101014500     145   CT
## 2       42      101  004202 1400000US42101004202 42101004202   42.02   CT
## 3       42      017  101805 1400000US42017101805 42017101805 1018.05   CT
## 4       42      091  200704 1400000US42091200704 42091200704 2007.04   CT
## 5       42      045  403001 1400000US42045403001 42045403001 4030.01   CT
## 6       42      101  004102 1400000US42101004102 42101004102   41.02   CT
## 7       42      101  000804 1400000US42101000804 42101000804    8.04   CT
## 8       42      101  025300 1400000US42101025300 42101025300     253   CT
## 9       42      101  028100 1400000US42101028100 42101028100     281   CT
## 10      42      101  008602 1400000US42101008602 42101008602   86.02   CT
##      ALAND AWATER                                                NAME.y
## 1   321465      0   Census Tract 145, Philadelphia County, Pennsylvania
## 2   453345      0 Census Tract 42.02, Philadelphia County, Pennsylvania
## 3  2641263      0      Census Tract 1018.05, Bucks County, Pennsylvania
## 4  4004196      0 Census Tract 2007.04, Montgomery County, Pennsylvania
## 5   492665      0   Census Tract 4030.01, Delaware County, Pennsylvania
## 6   391854      0 Census Tract 41.02, Philadelphia County, Pennsylvania
## 7   145445      0  Census Tract 8.04, Philadelphia County, Pennsylvania
## 8   581298      0   Census Tract 253, Philadelphia County, Pennsylvania
## 9   503541      0   Census Tract 281, Philadelphia County, Pennsylvania
## 10  369142      0 Census Tract 86.02, Philadelphia County, Pennsylvania
##      variable estimate moe rank                       geometry
## 1  B03002_001     1868 253    1 MULTIPOLYGON (((-75.15131 3...
## 2  B03002_001     5651 563    1 MULTIPOLYGON (((-75.15614 3...
## 3  B03002_001     2172 106    1 MULTIPOLYGON (((-75.13101 4...
## 4  B03002_001     3303 224    1 MULTIPOLYGON (((-75.30484 4...
## 5  B03002_001     3071 271    1 MULTIPOLYGON (((-75.30674 3...
## 6  B03002_001     8671 845    1 MULTIPOLYGON (((-75.16396 3...
## 7  B03002_001     4899 454    1 MULTIPOLYGON (((-75.17067 3...
## 8  B03002_001     4636 523    1 MULTIPOLYGON (((-75.18653 4...
## 9  B03002_001     4498 696    1 MULTIPOLYGON (((-75.15237 4...
## 10 B03002_001     2940 375    1 MULTIPOLYGON (((-75.22156 3...
popup <- paste0("GEOID: ", race_merged$GEOID, "<br>", "Number of Hispanic or Latino People: ", race_merged$estimate)
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = race_merged$estimate
)

map3<-leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = race_merged, 
              fillColor = ~pal(estimate), 
              color = "#b2aeae", # you need to use hex colors
              fillOpacity = 0.7, 
              weight = 1, 
              smoothFactor = 0.2,
              popup = popup) %>%
  addLegend(pal = pal, 
            values = race_merged$estimate, 
            position = "bottomright", 
            title = "Number of Hispanic or Latino People") 
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'
map3
readRenviron("~/.Renviron")
options(tigris_use_cache = TRUE)

# Get land area for tracts to calculate population density
sepabg <- tracts(state = "42", county = c("101", "017", "029", "045",
                                                "091"))
st_geometry(sepabg) <- NULL
head(sepabg)
##    STATEFP COUNTYFP TRACTCE       GEOID    NAME             NAMELSAD MTFCC
## 14      42      017  100302 42017100302 1003.02 Census Tract 1003.02 G5020
## 15      42      017  100304 42017100304 1003.04 Census Tract 1003.04 G5020
## 16      42      017  100306 42017100306 1003.06 Census Tract 1003.06 G5020
## 17      42      017  100307 42017100307 1003.07 Census Tract 1003.07 G5020
## 18      42      017  100401 42017100401 1004.01 Census Tract 1004.01 G5020
## 19      42      017  100402 42017100402 1004.02 Census Tract 1004.02 G5020
##    FUNCSTAT   ALAND  AWATER    INTPTLAT     INTPTLON
## 14        S 6861295 2033655 +40.0847565 -074.8915551
## 15        S 4018834   99694 +40.1122065 -074.8959059
## 16        S 1571083       0 +40.1061893 -074.8836008
## 17        S 2314482   71823 +40.0984599 -074.8984784
## 18        S 3245952    6352 +40.1348354 -074.8731389
## 19        S 3412871   20737 +40.1380006 -074.8558212
sepabg$A_KM2 <- sepabg$ALAND*1e-6

sf1_raw <- get_acs(geography = "tract", year = 2019, 
                         variables = c("B02001_001",  # Race total
                                       "B02001_002",   # Total Non-Hispanic White
                                       "B02001_003",   # Total Black population
                                       "B02001_004",   # Total Native American population 
                                       "B02001_005",   # Total Asian
                                       "B03002_001",  # Total Ethnicity
                                       "B03002_012"   # Total Hispanic or Latinx
                                       ),  
                         state = "PA", county = c("Philadelphia", "Bucks", 
                                                  "Chester", "Delaware", 
                                                  "Montgomery"))
## Getting data from the 2015-2019 5-year ACS
sf1 <- sf1_raw %>%
  group_by(GEOID) %>%
  summarise(race_tot = estimate[variable == "B02001_001"],
            white_count = estimate[variable == "B02001_002"],
            black_count = estimate[variable == "B02001_003"],
            asian_count = estimate[variable == "B02001_004"],
            native_count = estimate[variable == "B02001_005"],
            hispa_tot = estimate[variable == "B03002_001"],
            hispa_count = estimate[variable == "B03002_012"]
            ) %>%
  mutate(white_per = white_count/race_tot*100,
         black_per = black_count/race_tot*100,
         asian_per = asian_count/race_tot*100,
         native_per = native_count/race_tot*100,
         hispa_per = hispa_count/hispa_tot*100,
         pop_count = race_tot) %>%
  dplyr::select(GEOID, pop_count, white_per, black_per, asian_per, native_per, hispa_per)
knitr::kable(head(sf1), "simple")
GEOID pop_count white_per black_per asian_per native_per hispa_per
42017100102 2809 84.54966 2.064792 0.0000000 11.249555 2.634390
42017100103 2421 84.59314 4.337051 0.0413052 1.321768 13.878563
42017100104 4730 41.98732 21.754757 1.1205074 17.441861 19.978858
42017100105 3215 88.67807 1.648523 0.0000000 6.780715 8.087092
42017100201 4719 73.89277 6.166561 0.0000000 13.180759 5.340114
42017100206 5252 72.37243 7.159178 0.0000000 14.413557 7.806550
append_geoid <- function(address, geoid_type = 'tract') {

  if ("lat" %in% colnames(address) && "lon" %in% colnames(address)) {
    # Call for each row of the data
    geoids <- vector(mode="character", length = nrow(address))
    for (i in 1:nrow(address)) {
      geoids[i] <- call_geolocator_latlon(address$lat[i], address$lon[i])
    }
  } else {
    # If street, city, or state columns are factors, convert them
    # Call for each row of the data
    geoids <- vector(mode="character", length = nrow(address))
    for (i in 1:nrow(address)) {
      geoids[i] <- call_geolocator(
        as.character(address$street[i]),
        as.character(address$city[i]),
        as.character(address$state[i])
      )
    }
  }
  geoids
  # Append onto database
  address <- dplyr::mutate(address, geoid = geoids)

  # AABBBCCCCCCDEEE
  if (geoid_type == 'county') {
    end <- 5
  } else if (geoid_type == 'tract') {
    end <- 11
  } else if (geoid_type == 'block group') {
    end <- 12
  } else {
    end <- 15
  }
  address <- dplyr::mutate(address,
                           geoid = ifelse(is.na(.data$geoid), NA_character_, substr(.data$geoid, 1, end)))

  return(address)
}
final <- merge(sf1, census_final, by = "GEOID")
final[] <- lapply(final, function(x) {
    if(is.factor(x)) as.numeric(as.character(x)) else x
})
sapply(final, class)
##       GEOID   pop_count   white_per   black_per   asian_per  native_per 
## "character"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
##   hispa_per        pm25        pm10         no2          o3          co 
##   "numeric"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
##         so2          pb 
##   "numeric"   "numeric"
final <- final %>% dplyr::select(-GEOID, -pm10)
x <- cor(final, use = "complete.obs")
col <- colorRampPalette(c("darkorange", "white", "steelblue"))(20)
corrplot(x, type = "upper", order = "hclust", col = col)

final
##     pop_count  white_per   black_per  asian_per  native_per   hispa_per
## 1        2421 84.5931433  4.33705081 0.04130525  1.32176786 13.87856258
## 2        4730 41.9873150 21.75475687 1.12050740 17.44186047 19.97885835
## 3        5283 88.3399584  4.33465834 0.17035775  6.28430816  1.81714935
## 4        2555 80.4305284  5.28375734 0.31311155 12.91585127  2.66144814
## 5        2860 22.9370629  4.30069930 0.27972028 65.20979021  6.88811189
## 6        3385 68.4490399  7.09010340 0.00000000 16.33677991  5.31757755
## 7        2516 65.2225755 11.36724960 0.00000000 18.79968203  3.81558029
## 8        2543 39.3629571 32.91388124 2.24144711 14.31380260 11.01061738
## 9        1222 80.9328969  8.83797054 0.00000000  9.81996727  7.03764321
## 10       1514 88.3091149  0.72655218 0.00000000  8.45442536  6.07661823
## 11       2774 81.4708003  3.60490267 0.00000000  9.91348234  4.97476568
## 12       2056 76.0214008  4.47470817 0.00000000 16.48832685  2.91828794
## 13       2818 71.2207239  6.24556423 1.27750177 15.54293825  4.29382541
## 14       2670 91.9475655  2.50936330 0.00000000  3.10861423  2.09737828
## 15       3444 90.0696864  2.84552846 0.87108014  4.73286876  3.91986063
## 16       3750 73.9466667 10.40000000 0.00000000  9.94666667  7.01333333
## 17       2687 72.9065873  9.56457015 0.00000000 13.91886863  1.56308150
## 18       4038 78.6032689  6.33977216 0.32194156  7.70183259  6.09212481
## 19       4199 65.4917838 23.19599905 0.00000000  7.21600381  2.09573708
## 20       3227 71.5525256 13.91385187 0.00000000  9.91633096  6.91044314
## 21       2371 87.7266976  3.20539857 0.00000000  4.93462674  6.70603121
## 22       3114 84.8747592  7.70712909 0.00000000  0.57803468  4.91329480
## 23       3138 75.3027406 17.43148502 0.00000000  4.97131931  2.00764818
## 24       4039 62.4411983 31.31963357 0.00000000  4.53082446  6.21440951
## 25       1828 32.1663020 57.11159737 0.60175055  3.44638950  1.75054705
## 26       2221 38.2260243 55.29040973 0.00000000  3.06168393  4.27735254
## 27       2263 49.3592576 38.09102961 0.00000000  8.04242156  6.00972161
## 28       2918 75.4626456  2.91295408 0.00000000 18.88279644  7.64222070
## 29       4647 75.5756402  8.02668388 0.00000000 12.86851732  9.31783947
## 30       4243 58.0014141 18.26537827 0.00000000 16.63917040  8.17817582
## 31       4234 54.3221540 18.28058573 2.71610770 18.25696741 24.09069438
## 32       4419 85.8565286  2.73817606 0.00000000  3.96017198  8.10138040
## 33       4595 39.5647443  6.20239391 0.58759521 30.81610446 32.16539717
## 34       5596 70.1751251  6.95139385 0.00000000 15.29664046 15.35025018
## 35       3945 80.6844106  1.39416984 0.00000000  8.89733840 16.09632446
## 36       3959 43.8494569 19.29780248 0.00000000 32.78605709  8.86587522
## 37       2820 36.0283688 57.02127660 0.46099291  2.37588652  1.87943262
## 38       4426 26.9995481 60.09941256 1.46859467  3.36647085 10.32535020
## 39       4180 11.6985646 79.95215311 0.38277512  2.60765550 10.98086124
## 40       6203 44.3172658 41.83459616 0.00000000  6.93212961 15.75044333
## 41       7142 17.1800616 61.50938113 1.28815458 17.05404649  3.31839821
## 42       5859 17.4091142 65.28417819 2.35535074  7.35620413  8.70455709
## 43       4221 23.3120114 40.67756456 1.72944800 29.96920161  6.87040986
## 44       4015 88.0946451  0.49813200 0.00000000  9.19053549  7.82067248
## 45       6186 72.3892661  8.30908503 0.00000000 18.49337213  5.70643388
## 46       5800 96.3620690  0.29310345 0.41379310  1.70689655  1.34482759
## 47       4165 77.3829532  3.31332533 0.00000000  5.37815126 20.43217287
## 48       4757 85.7893630  4.07820055 0.00000000  7.50472987  2.87996637
## 49       5637 38.4424339 15.02572290 0.00000000 40.00354799 21.69593756
## 50       5326 83.5523845  0.88246339 0.00000000 14.73901615  3.88659407
## 51       5651 71.5802513  3.85772430 0.08847992 18.65156609  9.16651920
## 52       1577 16.6772353 70.32339886 0.00000000  8.18008878  7.67279645
## 53       6321 16.7062174 71.20708749 0.36386648  2.54706534 11.29568106
## 54       1158  5.6994819 88.16925734 0.00000000  0.34542314  3.02245250
## 55       6467 10.8705737 73.85186331 0.00000000  8.55110561  1.17519715
## 56       2759 10.2210946 81.76875680 0.00000000  7.68394346  0.36245016
## 57       4994  4.9058871 89.86784141 0.00000000  4.70564678  3.22386864
## 58       4512  3.0806738 81.16134752 0.08865248 11.45833333  5.98404255
## 59       4430  7.0428894 84.17607223 0.00000000  5.41760722  2.09932280
## 60       4426  1.2200633 97.62765477 0.00000000  0.00000000  0.36150023
## 61       3471  7.5482570 80.63958513 0.00000000  9.88187842  4.89772400
## 62       7082  7.0036713 81.12115222 3.62891838  6.91894945  0.50833098
## 63       3529  1.7852083 96.42958345 0.00000000  1.21847549  0.00000000
## 64       2595  3.5838150 92.98651252 0.00000000  1.00192678  2.42774566
## 65       4570  5.0109409 88.81838074 2.07877462  0.00000000  6.98030635
## 66       4466  3.1795790 87.48320645 0.00000000  4.63502015  1.47783251
## 67       2959 12.5042244 82.59547144 0.50692802  1.89253126  1.82494086
## 68       4720 10.5084746 80.36016949 0.00000000  1.27118644  1.86440678
## 69       1872 48.1837607 26.92307692 0.00000000 17.46794872  6.99786325
## 70       4163 62.1186644 23.73288494 0.00000000  7.20634158  6.19745376
## 71       4463 55.0302487 27.08940175 0.00000000  9.50033610  6.96840690
## 72       3807 25.1904387 66.79800368 0.00000000  2.44286840  4.88573680
## 73       3179  3.6174898 94.02327776 0.00000000  0.00000000  0.18873860
## 74       5354  8.0127008 91.42697049 0.00000000  0.00000000  3.67949197
## 75       5882  3.8252295 93.06358382 0.00000000  1.27507650  3.74022441
## 76       4628  1.9662921 94.87899741 0.17286085  1.14520311  1.53414002
## 77       5133  1.7533606 94.01909215 0.56497175  2.33781414  2.31833236
## 78       5050  2.0000000 93.92079208 0.21782178  1.20792079  1.76237624
## 79       6803  7.0263119 86.50595326 2.41070116  2.88108188  3.82184330
## 80       3312 51.2983092 27.29468599 0.87560386 14.16062802 10.92995169
## 81       2940 21.5986395 66.49659864 0.37414966  8.43537415  3.87755102
## 82       3667 55.6858467 18.24379602 1.52713390 20.34360513  6.05399509
## 83       2622 54.9961861 10.45003814 0.00000000 25.51487414  9.07704043
## 84       2289 51.9440804 11.88291830 0.87374399 23.46002621 10.96548711
## 85       6303 50.8488022 21.94193241 0.00000000 22.57655085  9.18610186
## 86       6201 68.5857120  7.83744557 0.19351717 19.60974036  7.61167554
## 87       2958 35.2603110 36.24070318 0.91277890 22.44759973  3.71872887
## 88       3345 17.6382661 65.23168909 0.20926756  9.08819133  3.91629297
## 89       4721  5.8885829 87.60855751 0.00000000  3.28320271  2.96547342
## 90       4252  1.8814675 96.35465663 0.00000000  0.00000000  2.30479774
## 91       3261  5.5504446 93.49892671 0.00000000  0.00000000  3.80251457
## 92       4477  0.2903730 94.21487603 0.49140049  0.17869109 10.58744695
## 93       2318 13.1147541 82.83002588 0.00000000  1.72562554  5.43572045
## 94       6190 10.3231018 86.28432956 0.00000000  0.56542811  2.40710824
## 95       4863  5.8400165 89.01912400 0.00000000  2.30310508  1.89183632
## 96       6925  1.7617329 96.33212996 0.30324910  0.00000000  0.30324910
## 97       3344  3.5885167 94.94617225 0.00000000  0.00000000  5.32296651
## 98       2816  1.5269886 96.09375000 0.00000000  1.49147727  1.13636364
## 99       3827  5.2260256 91.06349621 0.00000000  0.00000000  1.56780768
## 100      4544  6.3600352 89.37059859 0.77024648  1.36443662  2.81690141
## 101      1704  4.6361502 86.73708920 2.99295775  3.16901408  2.99295775
## 102      3612  5.4540421 92.74640089 0.00000000  0.16611296  1.93798450
## 103      4258  7.9614843 84.78158760 1.43259746  3.75763269  1.87881635
## 104      2790 18.1362007 72.86738351 0.50179211  5.16129032  5.16129032
## 105      3549  5.0154973 81.09326571 0.22541561  2.19780220  2.90222598
## 106      3704  4.4816415 91.14470842 0.00000000  0.05399568  1.05291577
## 107      5939  3.1486782 93.11331874 0.00000000  0.00000000  6.51624853
## 108      3344  1.8540670 98.14593301 0.00000000  0.00000000  1.58492823
## 109      6864  1.3840326 97.08624709 0.00000000  0.23310023  1.15093240
## 110      4472  2.0796064 95.99731664 0.67084079  0.00000000  1.47584973
## 111      1957 44.8645887 38.98824732 0.00000000 14.97189576  2.75932550
## 112      5618  2.7589890 92.93342827 0.00000000  1.15699537  3.09718761
## 113      5161  2.7126526 94.03216431 0.00000000  0.00000000  0.73629142
## 114      1890 31.0582011 58.14814815 0.00000000  8.57142857  4.28571429
## 115      3189 19.4731891 74.31796802 0.00000000  3.51207275  3.76293509
## 116      2776 33.2853026 45.49711816 3.96253602 13.54466859  3.45821326
## 117       633 23.6966825 24.17061611 0.00000000 38.54660348  9.79462875
## 118      4265 15.7561547 75.26377491 1.10199297  4.47831184  2.62602579
## 119      1765 23.9093484 61.98300283 0.00000000  7.53541076  4.47592068
## 120      2778 26.3858891 59.53923686 1.65586753  7.37940965  9.86321094
## 121      3100 57.6129032 26.90322581 0.00000000 10.87096774  3.70967742
## 122      2504 88.5383387  5.99041534 0.00000000  4.23322684  0.00000000
## 123      3088 82.8044041  6.47668394 0.71243523  4.82512953 11.94948187
## 124      3987 68.7735139 25.30724856 0.00000000  3.78730875  6.79709054
## 125      2832 83.2980226  3.67231638 0.00000000 10.16949153  2.68361582
## 126      3827 78.2336033  5.98379932 0.00000000  6.68931278  0.88842435
## 127      1972 16.5821501 70.38539554 0.00000000  1.82555781  8.36713996
## 128      1757 29.3682413 62.43597040 0.51223677  0.51223677  7.68355151
## 129      3888 24.3827160 61.26543210 0.00000000  7.71604938  5.96707819
## 130      2659 21.5494547 66.52877021 0.48890560  7.33358405  5.41556976
## 131      1935 86.9250646  3.61757106 0.00000000  6.14987080  5.37467700
## 132      4470 64.4071588 12.32662192 0.29082774  4.98881432 30.80536913
## 133      1868 11.4025696 73.28693790 0.00000000  1.44539615 22.48394004
## 134      4305 35.5865273 47.03832753 0.00000000 11.08013937  5.90011614
## 135      2686 37.7513031 49.73938943 0.00000000  7.40878630  8.67460908
## 136       804  2.2388060 82.33830846 0.00000000 12.68656716  1.24378109
## 137      3891  6.7077872 84.24569519 0.00000000  0.77101002 12.72166538
## 138      2279  4.3878894 88.45985081 5.74813515  0.35103115  7.72268539
## 139      4242  0.1650165 95.56812824 0.09429514  2.09806695  4.24328147
## 140      5347  3.6843090 93.92182532 0.26182906  1.68318683  4.76902936
## 141      4086 50.7831620 33.70044053 1.44395497 10.67058248  5.65345081
## 142      1596 45.1127820 26.37844612 1.87969925  3.13283208 67.04260652
## 143      2587 57.6729803 11.63509857 0.19327406  1.50753769 34.20950908
## 144      6969 92.0505094  2.84115368 0.00000000  3.07074186  4.57741426
## 145      6203 63.7594712  6.14218926 1.12848622 11.97807512 23.92390779
## 146      2098 40.7530982 22.83126787 0.00000000  7.62631077 69.68541468
## 147      3005 48.5191348 12.71214642 0.00000000  5.65723794 74.57570715
## 148      4633 26.8292683 33.75782430 1.16555148  0.00000000 58.21282107
## 149      3180  7.4213836 89.33962264 0.15723270  1.54088050 10.59748428
## 150      2241 15.6180277 76.30522088 0.00000000  1.33868809  9.72780009
## 151      2809  1.7443930 94.80242079 0.00000000  0.00000000  1.85119260
## 152      2704  7.9881657 91.71597633 0.00000000  0.29585799  2.95857988
## 153      2965  1.1804384 96.35750422 0.00000000  0.00000000  0.43844857
## 154      3092  1.5200517 98.47994825 0.00000000  0.00000000  1.22897801
## 155      4412  0.7252947 97.91477788 0.00000000  0.00000000  0.88395286
## 156      2884 23.1276006 68.16920943 0.52011096  1.66435506  3.81414702
## 157      3565  1.1781206 91.89340813 0.00000000  0.00000000  7.01262272
## 158      2835  4.4797178 75.69664903 0.00000000  7.12522046  2.32804233
## 159      3842  0.2342530 95.96564289 0.80687142  2.99323269  0.00000000
## 160      2538  5.5555556 91.01654846 0.74862096  0.03940110  5.47675335
## 161      2392 11.7056856 80.85284281 0.00000000  0.54347826 14.84113712
## 162      6943 33.2421144 18.86792453 0.87858275  1.02261270 77.74737145
## 163      4972 31.7980692 13.21399839 0.48270314  0.08045052 88.85760257
## 164      4603 41.9509016 13.53465131 0.26069954  0.56484901 83.55420378
## 165      3203 36.2472682 31.59537933 0.62441461  0.31220731 55.85388698
## 166      5188 29.7224364 22.55204318 1.04086353  0.00000000 78.48882035
## 167      5862 46.4346639 29.68270215 2.37120437  0.71647902 44.62640737
## 168      6279 43.2075171 19.22280618 0.00000000  0.35037426 56.63322185
## 169      3253 86.1973563  1.56778358 0.00000000  1.96741469 11.77374731
## 170      6038 91.4706857  0.01656178 0.00000000  4.42199404  8.52931434
## 171      4005 97.2784020  0.49937578 0.00000000  0.00000000  1.92259675
## 172      2242 98.5280999  0.04460303 0.22301517  0.93666369  5.79839429
## 173      7750 28.5290323 16.92903226 0.00000000  2.82580645 69.22580645
## 174      8194 42.4334879 15.34049304 0.48816207 13.07053942 63.27800830
## 175      8721 42.0479303 21.06409815 0.00000000  4.01330123 65.14161220
## 176      4432 43.7951264  7.85198556 2.18862816  0.00000000 88.06407942
## 177      3732 48.3386924  8.97642015 0.00000000  0.00000000 80.49303323
## 178      7064 26.1041903 40.81257078 0.35390713  0.00000000 58.67780294
## 179      6878 16.6763594 36.21692352 1.77377145  5.33585345 61.35504507
## 180      4498 28.8572699 39.21742997 0.00000000  0.00000000 64.25077812
## 181      1320  6.5151515 72.87878788 0.00000000  1.36363636 19.62121212
## 182      3552  5.6306306 89.07657658 0.33783784  2.05518018  5.71509009
## 183      3415  1.0834553 98.30161054 0.38067350  0.00000000  4.24597365
## 184      4804  6.5986678 88.92589509 0.00000000  0.16652789  5.97418818
## 185      2842  2.3926812 81.70302604 0.91484870  0.35186488 14.07459536
## 186      2744  1.6399417 95.59037901 0.00000000  0.00000000  3.24344023
## 187      2578  3.0256012 91.15593483 0.00000000  4.80993018  2.36617533
## 188      1695 66.4306785 19.23303835 0.00000000  6.84365782 10.08849558
## 189      2289 49.9344692 34.33813893 0.00000000 12.66928790  3.84447357
## 190      3085 83.0794165 11.28038898 1.19935170  2.39870340  4.70016207
## 191      5342 84.8558592 10.27704979 0.50542868  1.96555597  4.41782104
## 192      2481 85.8524788  8.22249093 0.00000000  2.29746070  3.78879484
## 193      2419 72.7159983 19.80157090 0.00000000  2.56304258  7.60644895
## 194      3485 84.5050215  9.03873745 0.71736011  2.15208034  5.68149211
## 195      3980 83.3919598 10.37688442 0.37688442  1.08040201  5.20100503
## 196      3879 80.6651199 11.75560712 0.51559680  3.22248002  2.88734210
## 197      2372 70.6998314 17.70657673 0.00000000  5.26981450  3.45699831
## 198      5727 77.7894185 12.08311507 0.00000000  3.37000175  3.49222979
## 199      5141 56.1563898 31.31686442 0.00000000  4.80451274  3.46236141
## 200      1434 93.3054393  3.48675035 0.00000000  1.67364017  2.85913529
## 201      1534 88.7222947  8.53976532 0.00000000  2.73794003  3.91134289
## 202      1275 94.1960784  0.86274510 0.62745098  0.70588235  4.86274510
## 203      1300 66.6153846 27.76923077 0.00000000  1.61538462  0.15384615
## 204      2391 60.2258469 32.70598076 0.00000000  1.46382267  2.04935174
## 205      4463 29.8678019 59.01859736 0.00000000  1.50123235  4.16760027
## 206      4633 27.8653141 62.18432981 0.00000000  0.84178718  5.71983596
## 207      1634 31.4565483 53.73317013 0.00000000 12.79069767  6.67074663
## 208      3659 36.1847499 55.09702104 0.00000000  3.58021317  5.35665482
## 209      1445 19.9307958 72.24913495 0.00000000  0.62283737  3.59861592
## 210      3807 14.2369320 83.18886262 0.00000000  1.33963751  2.83687943
## 211      3895 13.0680359 81.43774069 0.56482670  1.95121951  0.82156611
## 212      3230 10.5263158 83.68421053 1.82662539  0.00000000 11.54798762
## 213      3823  8.4488622 90.00784724 0.00000000  0.73240910  0.62777923
## 214      2684 11.1028316 83.56929955 0.00000000  1.22950820  3.05514158
## 215      4293  6.2893082 90.58933147 0.20964361  0.32611228  1.81691125
## 216      2232  2.5537634 92.29390681 0.00000000  0.00000000  3.00179211
## 217      2737  1.5710632 96.96748265 0.07307271  0.00000000  0.18268177
## 218      7793  5.7230848 91.37687668 1.15488259  0.80841781  4.28589760
## 219      4636 11.4754098 84.18895600 0.00000000  2.78257118  0.04314064
## 220      4222 12.4111795 79.67787778 0.00000000  0.80530554  2.67645666
## 221      2764 38.7843705 48.26338640 0.00000000  3.40086831  2.74963821
## 222      2795 53.4168157 36.77996422 0.00000000  0.96601073  3.22003578
## 223      3874 48.7351575 35.85441404 0.00000000  6.50490449  9.03458957
## 224      1778  5.7367829 91.78852643 0.00000000  1.40607424  0.00000000
## 225      4409  1.9505557 95.07824904 0.00000000  0.27217056  0.52166024
## 226      3055  3.9279869 90.54009820 0.52373159  0.32733224  2.48772504
## 227      3192  3.2268170 95.39473684 0.00000000  0.00000000  0.93984962
## 228      3791  3.4819309 95.04088631 0.21102611  0.18464785  1.76734371
## 229      3488  1.4048165 95.61353211 0.00000000  0.02866972  1.69151376
## 230      4927  3.0850416 95.90014207 0.00000000  1.01481632  0.87274203
## 231      5660  1.0954064 98.51590106 0.00000000  0.00000000  0.45936396
## 232      5065  1.6979269 94.43237907 0.00000000  0.00000000  1.53998026
## 233      6528  0.2297794 95.17463235 0.00000000  0.00000000  2.51225490
## 234      6987  0.5868041 96.04980678 0.00000000  0.00000000  0.65836554
## 235      4189  8.6655526 77.96610169 0.00000000  3.17498210 11.12437336
## 236      2147 21.5183978 65.20726595 0.00000000  0.32603633 11.50442478
## 237      2635 10.3605313 79.73434535 0.00000000  6.60341556  0.64516129
## 238      2806  6.7355666 62.33071989 0.00000000 20.34925160  9.19458304
## 239      4647 10.3077254 76.35033355 0.00000000  6.45577792  8.43554982
## 240      6019  6.9612893 66.80511713 0.44857950 17.62751288  4.90114637
## 241      2975 11.0252101 56.80672269 0.00000000 20.73949580 17.44537815
## 242      6825  8.2051282 56.73260073 0.00000000  4.82051282 34.37362637
## 243      5133 10.1694915 75.06331580 0.50652640  7.71478667 15.05941944
## 244      4163  0.7446553 94.30699015 0.00000000  1.51333173  2.52221955
## 245      4960  1.5322581 88.46774194 0.00000000  2.96370968  0.48387097
## 246      4564  3.2427695 94.82909728 0.15337423  0.83260298  2.05959684
## 247      4201  1.9519162 92.04951202 0.00000000  0.49988098  1.42823137
## 248      4075 38.7239264 52.58895706 0.00000000  4.71165644  3.65644172
## 249      4111  1.1675991 96.13232790 0.00000000  1.80004865  0.87569934
## 250      4498  7.5811472 85.90484660 0.00000000  0.86705202  9.29301912
## 251      5184 17.8433642 69.17438272 0.00000000 11.22685185 17.09104938
## 252      6738  0.1038884 77.20391808 0.00000000 13.60937964  1.05372514
## 253      3873  9.1660212 82.49419055 0.54221534  5.88690937 10.53446940
## 254      2480 14.3145161 51.08870968 0.00000000 14.39516129 33.62903226
## 255      6583 10.0106334 68.37308218 0.00000000  7.01807686 21.28209023
## 256      2577 22.6232053 35.04074505 0.69848661  2.71633683 61.46682189
## 257      4619 23.1651873 17.75276034 0.00000000 14.22385798 62.15631089
## 258      4056 30.3254438 21.96745562 5.84319527 10.37968442 69.45266272
## 259      7043 23.4985092 30.24279426 0.97969615  9.34260968 53.84069289
## 260      8124 20.7779419 37.96159527 0.00000000  9.17035943 46.17183653
## 261      5158 12.8538193 38.30942226 0.27142303 26.28925940 36.25436216
## 262      4917 19.3817368 38.88549929 2.98962782  2.42017490 47.14256661
## 263      3170 28.5173502 23.75394322 2.39747634  0.00000000 57.60252366
## 264      3874 24.5482705 51.91017037 0.82601962  8.98296335 30.97573567
## 265      4830 38.6542443 42.07039337 0.00000000  2.85714286 36.64596273
## 266      4324 24.3987049 52.47456059 3.76965772  1.98889917 31.56799260
## 267      7629 19.8584349 62.60322454 0.24904968  3.18521431 23.39756193
## 268      5984 30.2473262 47.76069519 0.08355615  1.38703209 29.66243316
## 269      8912 17.5044883 58.23608618 0.76301616  2.15439856 30.30745063
## 270      6079 12.4691561 61.52327685 1.08570489  5.11597302 23.86905741
## 271      7254 25.9718776 55.27984560 0.00000000  6.43782741 26.04080507
## 272      6988 38.2083572 35.90440756 0.38637665  9.55924442 23.48311391
## 273      4238 47.3100519 38.83907504 0.00000000  6.81925437 17.41387447
## 274      5570 45.0987433 26.10412926 0.25134650 15.06283662 15.45780969
## 275      4553 33.2747639 41.81858116 0.00000000 12.12387437 23.54491544
## 276      7694 41.7598128 24.04471016 4.43202495 18.22199116 22.92695607
## 277      5274 24.3458476 45.43041335 0.00000000  9.65111870 25.19908987
## 278      4346 15.7156006 53.70455591 2.57708237 10.07823286 34.49148642
## 279      5032 27.5437202 49.86089030 0.23847377  8.98251192 24.42368839
## 280      7632 40.8936059 20.70230608 0.64203354 20.71540881 25.51100629
## 281      7061 31.7093896 21.20096304 0.00000000 21.55502054 23.77850163
## 282      5921 36.1763216 13.10589427 0.28711366 34.06519169 29.21803749
## 283      6643 67.9211200  4.84720759 0.00000000 22.97154900 11.78684329
## 284      4337 50.2421028 20.49804012 0.00000000 17.54669126 12.84297902
## 285      5998 64.1547182 16.18872958 0.00000000 10.17005669 18.97299100
## 286      5175 46.4347826 28.11594203 0.00000000  5.85507246 24.40579710
## 287      4056 29.9063116 43.86094675 0.00000000  3.67357002 30.22682446
## 288      5641 39.3724517 34.03651835 3.10228683  4.23683744 30.93423152
## 289      7382 67.9355188 15.52424817 0.00000000  3.67109185 27.64833378
## 290      4070 52.4815725 20.04914005 0.00000000  3.85749386 34.05405405
## 291      3619 52.3901630 36.66758773 0.00000000  0.00000000 25.67007461
## 292      6143 53.5406153 20.02279017 0.00000000 10.76021488 18.19957675
## 293      7611 67.0345552  8.89502037 0.00000000 10.48482460 16.04256996
## 294      3701 61.3617941 18.69764928 0.00000000  3.64766279 23.88543637
## 295      8394 58.8277341 21.86085299 0.20252561  4.13390517 28.30593281
## 296      5064 83.3925750  7.54344392 0.00000000  2.48815166 11.23617694
## 297      4048 85.5978261  7.80632411 0.00000000  1.48221344 10.82015810
## 298      2960 64.4594595  6.62162162 0.00000000 18.98648649 18.75000000
## 299      4147 75.1627683 13.38316856 0.00000000  5.90788522 11.45406318
## 300      5619 49.7063534 21.71204841 0.00000000 24.82648158  6.05089874
## 301      3975 47.2201258 15.62264151 0.00000000 20.42767296 18.99371069
## 302      7288 67.3161361 20.82875960 0.00000000  7.08013172 13.78979144
## 303      6589 63.3935347 19.10760358 0.00000000  7.46699044 21.65730763
## 304      4860 58.0452675  9.19753086 0.00000000 25.32921811  5.88477366
## 305      6380 70.1097179 10.21943574 0.32915361  2.64890282 19.74921630
## 306      3133 63.8365784 16.53367380 0.00000000 11.77784871 15.22502394
## 307      3082 75.5678131 11.19402985 0.00000000  9.86372485 16.77482154
## 308      5729 79.0713912  9.39081864 0.00000000  9.93192529 16.58230058
## 309      3354 76.8932618 18.42576029 0.00000000  3.19022063  5.21765057
## 310      8044 71.5440080  6.81253108 0.29835903 19.43063153  5.38289408
## 311      4243 55.9274099 17.29908084 0.00000000 15.93212350 17.20480792
## 312      5672 65.0564175 22.14386460 0.00000000  8.65655853 12.62341326
## 313      2496 48.3974359 28.96634615 0.24038462 13.42147436 20.07211538
## 314      6879 57.2612298 23.31734264 0.04361099 10.85913650  9.12923390
## 315      3835 69.4393742  8.47457627 0.00000000 10.37809648 12.46414602
## 316      4600 89.8260870  4.78260870 0.00000000  4.36956522  7.43478261
## 317      5485 59.5624430 29.29808569 0.51048314  0.00000000 11.19416591
## 318      4309 78.9278255  3.15618473 0.27848689  4.50220469 14.69018334
## 319      5316 56.5086531 20.20316027 0.00000000  4.81565087 23.79608728
## 320      4187 85.4549797  8.78910915 0.00000000  5.25435873  3.98853594
## 321      4686 92.7870252  2.68886044 0.00000000  3.13700384  4.26803244
## 322      5550 88.8468468  1.60360360 0.00000000  2.41441441 12.63063063
## 323      4516 85.3410097  9.38883968 0.22143490  2.10363153  6.64304694
## 324      7948 71.5525918  9.48666331 0.00000000 12.12883744  5.76245596
## 325      5272 79.4385432  2.65553869 0.17071320 16.82473445  2.40895296
## 326      3661 76.3725758  1.99399071 0.00000000 19.31166348  8.38568697
## 327      4931 56.1752180  6.10423849 0.00000000 27.86453052  5.51612249
## 328      4158 71.5728716  1.03415103 0.00000000 18.61471861  5.02645503
## 329      6806 77.4610638  2.24801646 0.51425213 12.66529533  2.17455187
## 330      5548 71.5032444  1.82047585 0.77505407 16.76279740 10.61643836
## 331      3472 56.2788018 11.49193548 0.23041475 27.30414747  3.71543779
## 332      4049 67.6957273 20.69646826 0.00000000  2.34625834 15.18893554
## 333      5153 88.2592664  5.97710072 0.00000000  0.00000000  8.13118572
## 334      6239 84.9174547  6.81198910 0.00000000  3.92691136 11.07549287
## 335      5111 91.4693798  3.03267462 1.64351399  1.85873606  8.04147916
## 336      3969 80.5744520  0.32753842 0.32753842 13.85739481  6.14764424
## 337      3705 72.6315789 22.56410256 0.00000000  0.00000000  9.47368421
## 338      6592 91.0497573  6.76577670 0.00000000  1.44114078  1.88106796
## 339      1037 78.9778206 12.82545805 0.57859209  1.44648023 11.28254581
## 340      5331 79.8724442 11.14237479 0.00000000  5.47739636  5.58994560
## 341      4213 88.3930691  0.16615239 0.00000000  5.38808450  7.57180157
## 342      2604 72.6958525  5.26113671 0.00000000 11.52073733 12.17357911
## 343      3612 81.7552602  8.36101883 0.24916944  4.98338870  8.72093023
## 344      4702 56.9757550 10.27222459 0.00000000 26.90344534  3.38153977
## 345      5759 81.3335649 10.14064942 1.02448342  7.32766105 12.95363778
## 346      3736 47.7516060 42.26445396 0.00000000  4.41648822  7.97644540
## 347      3287 46.8512321 18.92303012 0.60845756 26.01156069 10.92181320
## 348      5313 39.6009787 41.69019386 0.00000000  8.46979108 19.63109354
## 349      2474 90.1778496  5.73969281 0.48504446  1.45513339  6.83104285
## 350      5006 95.8250100  0.11985617 0.19976029  1.27846584  5.37355174
## 351      2431 71.9457014  9.46112711 1.27519539  1.23406006 38.66721514
## 352       845 64.8520710 21.06508876 0.00000000  0.35502959 22.36686391
## 353      3067 58.6566678 24.09520704 0.00000000  2.11933485 26.27975220
## 354      2484 85.4267311  9.70209340 0.00000000  2.85829308  1.81159420
## 355      1830 83.1693989  5.19125683 0.00000000  6.12021858  0.76502732
## 356      1395 80.5017921 13.33333333 0.50179211  4.08602151  2.93906810
## 357      2816 75.4261364 15.34090909 0.00000000  5.07812500  2.80539773
## 358      3636 63.0088009 27.11771177 0.44004400  2.50275028  3.52035204
## 359      3251  5.5982775 90.37219317 0.00000000  0.00000000  3.72193171
## 360        40 62.5000000  0.00000000 7.50000000  0.00000000  0.00000000
## 361       554 89.5306859  8.48375451 0.00000000  1.98555957  5.77617329
## 362         0        NaN         NaN        NaN         NaN         NaN
## 363         0        NaN         NaN        NaN         NaN         NaN
## 364         0        NaN         NaN        NaN         NaN         NaN
## 365         0        NaN         NaN        NaN         NaN         NaN
## 366         0        NaN         NaN        NaN         NaN         NaN
## 367      1811 14.7432358 74.15792380 2.04307013  0.88348978 19.32633904
##         pm25      no2         o3        co      so2          pb
## 1   7.876399 24.10012 0.03629215 0.5802992 1.163807 0.005398372
## 2   7.872433 24.03197 0.03621404 0.5752292 1.167688 0.005411446
## 3   7.895433 22.58762 0.03447842 0.5563563 1.243126 0.003981212
## 4   8.062386 22.56431 0.03421509 0.5980918 1.252842 0.002815069
## 5   7.634107 22.39948 0.03424346 0.5227076 1.305566 0.006580307
## 6   7.671309 22.47043 0.03424054 0.5529315 1.324246 0.006245931
## 7   7.662028 22.48747 0.03425637 0.5593391 1.337441 0.006616675
## 8   7.631549 22.40438 0.03425929 0.5231976 1.316991 0.006944825
## 9   7.626709 22.40415 0.03427090 0.5218728 1.323791 0.007214087
## 10  7.625414 22.49177 0.03427469 0.5596542 1.348397 0.007221687
## 11  7.618790 22.47242 0.03427302 0.5506557 1.345511 0.007345163
## 12  7.619603 22.42344 0.03427579 0.5289957 1.334370 0.007420908
## 13  7.621876 22.39748 0.03428038 0.5180747 1.326588 0.007401185
## 14  7.620745 22.37353 0.03429081 0.5076143 1.320767 0.007454978
## 15  7.624206 22.34797 0.03429703 0.4976318 1.310013 0.007371528
## 16  7.605289 22.42369 0.03428681 0.5277413 1.340564 0.007696623
## 17  7.612378 22.39871 0.03429144 0.5171353 1.332944 0.007642275
## 18  7.592708 22.48680 0.03428508 0.5561374 1.353303 0.007640434
## 19  7.564281 22.45583 0.03429573 0.5403780 1.354269 0.008066403
## 20  7.596746 22.40773 0.03430018 0.5195600 1.340629 0.007875011
## 21  7.613951 22.36216 0.03431374 0.5003345 1.324150 0.007714263
## 22  7.605507 22.36765 0.03432145 0.5010936 1.329482 0.007856120
## 23  7.584252 22.41067 0.03430657 0.5197655 1.344708 0.008022660
## 24  7.533403 22.45578 0.03430572 0.5389035 1.358809 0.008332543
## 25  7.478217 22.48517 0.03431647 0.5510116 1.366124 0.008621867
## 26  7.487866 22.46664 0.03431547 0.5422846 1.364546 0.008612622
## 27  7.507121 22.44787 0.03431721 0.5335197 1.361557 0.008542527
## 28  7.533780 22.42412 0.03432297 0.5224638 1.356104 0.008425964
## 29  7.563998 22.41164 0.03431788 0.5183133 1.349308 0.008218301
## 30  7.588802 22.37671 0.03433205 0.5025419 1.336908 0.008051579
## 31  7.547997 22.39745 0.03434706 0.5073237 1.349473 0.008370581
## 32  7.566162 22.38408 0.03435109 0.5015613 1.342977 0.008249101
## 33  7.535576 22.40634 0.03434312 0.5115221 1.353198 0.008443775
## 34  7.507603 22.42055 0.03434228 0.5173491 1.358949 0.008593828
## 35  7.492054 22.43252 0.03433560 0.5235954 1.362023 0.008662142
## 36  7.446027 22.45010 0.03433953 0.5304290 1.367726 0.008876944
## 37  7.473175 22.44881 0.03432866 0.5319334 1.365183 0.008740092
## 38  7.439483 22.46614 0.03433094 0.5393651 1.368769 0.008887066
## 39  7.430947 22.48343 0.03433089 0.5477032 1.369755 0.008907473
## 40  7.495382 22.50746 0.03432495 0.5612180 1.365518 0.008411709
## 41  7.406523 22.49122 0.03434931 0.5481125 1.371537 0.009026912
## 42  7.404788 22.46985 0.03434496 0.5384184 1.371503 0.009054836
## 43  7.393603 22.47174 0.03435365 0.5374928 1.372337 0.009105709
## 44  7.397886 22.47227 0.03436981 0.5339922 1.371919 0.009085820
## 45  7.423661 22.45270 0.03435424 0.5284888 1.369634 0.008978701
## 46  7.421454 22.45634 0.03437292 0.5257474 1.369584 0.008985561
## 47  7.461286 22.43571 0.03435537 0.5208599 1.365257 0.008818074
## 48  7.449365 22.44106 0.03436813 0.5201963 1.366360 0.008868300
## 49  7.501287 22.41772 0.03436010 0.5124289 1.358680 0.008630815
## 50  7.539218 22.39841 0.03436707 0.5034081 1.349441 0.008423838
## 51  7.530253 22.40619 0.03438021 0.5032168 1.350639 0.008467716
## 52  7.621657 22.53780 0.03446370 0.5420512 1.340469 0.007409541
## 53  7.605328 22.53299 0.03443179 0.5508802 1.348014 0.007514842
## 54  7.595177 22.53091 0.03444144 0.5459811 1.347767 0.007640311
## 55  7.579370 22.52845 0.03441707 0.5526356 1.353035 0.007746597
## 56  7.544567 22.52325 0.03440328 0.5533530 1.358147 0.008047880
## 57  7.570977 22.52789 0.03440006 0.5573150 1.355851 0.007797513
## 58  7.600982 22.53228 0.03440336 0.5593314 1.352346 0.007508764
## 59  7.623622 22.53526 0.03439890 0.5626377 1.350488 0.007270041
## 60  7.659705 22.53898 0.03437463 0.5711122 1.349675 0.006854431
## 61  7.588789 22.53097 0.03437890 0.5647250 1.356179 0.007595501
## 62  7.540987 22.52421 0.03438247 0.5594968 1.360201 0.008046938
## 63  7.582576 22.53017 0.03435839 0.5685828 1.358351 0.007634360
## 64  7.623483 22.53496 0.03435703 0.5719248 1.354981 0.007224582
## 65  7.661830 22.53904 0.03435592 0.5748352 1.351633 0.006820153
## 66  7.711981 22.54403 0.03435679 0.5780874 1.346653 0.006276217
## 67  7.682815 22.54066 0.03434207 0.5782462 1.350935 0.006582231
## 68  7.598425 22.53073 0.03433798 0.5725110 1.358118 0.007464476
## 69  7.640609 22.53285 0.03431740 0.5767480 1.355270 0.007016764
## 70  7.672758 22.53855 0.03432831 0.5788531 1.352616 0.006681496
## 71  7.733498 22.54458 0.03432339 0.5833505 1.347741 0.006014899
## 72  7.760607 22.54794 0.03433620 0.5839565 1.344264 0.005730847
## 73  7.753485 22.54782 0.03435474 0.5810947 1.342568 0.005818055
## 74  7.790104 22.55082 0.03434635 0.5845401 1.339788 0.005414857
## 75  7.772927 22.54984 0.03436852 0.5801004 1.337836 0.005606944
## 76  7.842658 22.55540 0.03435934 0.5858364 1.330329 0.004851182
## 77  7.837298 22.55481 0.03435157 0.5866335 1.332922 0.004911981
## 78  7.850593 22.55567 0.03433910 0.5890043 1.333827 0.004773957
## 79  7.835374 22.55408 0.03432479 0.5895718 1.337876 0.004927753
## 80  7.774344 22.54785 0.03431381 0.5865611 1.344568 0.005558214
## 81  7.823248 22.55204 0.03430693 0.5899692 1.340406 0.005029946
## 82  7.753079 22.54447 0.03430430 0.5853270 1.346323 0.005773444
## 83  7.748830 22.54274 0.03429752 0.5848999 1.346298 0.005803395
## 84  7.721535 22.53463 0.03428373 0.5813324 1.346518 0.006062434
## 85  7.719357 22.53752 0.03429309 0.5821741 1.348037 0.006116755
## 86  7.765567 22.53068 0.03425472 0.5812981 1.335028 0.005380278
## 87  7.802992 22.54262 0.03426723 0.5871790 1.337277 0.005084968
## 88  7.851096 22.55220 0.03428315 0.5919157 1.336907 0.004683547
## 89  7.906014 22.55951 0.03431082 0.5947054 1.331758 0.004220389
## 90  7.909127 22.56005 0.03432546 0.5937381 1.328906 0.004194407
## 91  7.924345 22.56126 0.03433331 0.5937692 1.324601 0.004051380
## 92  7.909682 22.56029 0.03434509 0.5915943 1.323552 0.004178203
## 93  7.914282 22.56192 0.03438257 0.5860423 1.307418 0.004074820
## 94  7.901384 22.56122 0.03438605 0.5848245 1.309658 0.004202887
## 95  7.929224 22.56184 0.03435370 0.5914471 1.316431 0.003978757
## 96  7.960702 22.56390 0.03433619 0.5953279 1.316632 0.003713956
## 97  7.959938 22.56379 0.03431746 0.5972010 1.323241 0.003736493
## 98  7.972793 22.56463 0.03430455 0.5988091 1.324161 0.003629227
## 99  7.948752 22.56248 0.03429523 0.5978891 1.328242 0.003834652
## 100 7.965982 22.56340 0.03428262 0.5991819 1.326476 0.003687124
## 101 7.924591 22.55894 0.03427375 0.5967198 1.329846 0.003976915
## 102 7.991236 22.56446 0.03426215 0.6007646 1.322132 0.003423465
## 103 7.917901 22.55544 0.03425150 0.5955331 1.325102 0.003936902
## 104 7.857522 22.54425 0.03423926 0.5895867 1.324003 0.004339650
## 105 8.055573 22.56969 0.03424893 0.6045472 1.312952 0.002962680
## 106 8.058413 22.57091 0.03428471 0.6043025 1.313399 0.002993469
## 107 8.002799 22.56687 0.03431515 0.5996174 1.315884 0.003382747
## 108 8.029546 22.56871 0.03431396 0.6009986 1.309367 0.003174749
## 109 8.000488 22.56657 0.03433325 0.5974134 1.306572 0.003373458
## 110 7.957224 22.56378 0.03435246 0.5929264 1.309690 0.003716183
## 111 8.042265 22.56908 0.03432008 0.5999244 1.293230 0.003049929
## 112 8.053202 22.57021 0.03431043 0.6021060 1.300797 0.002999439
## 113 8.086738 22.57275 0.03429215 0.6051856 1.298781 0.002796355
## 114 8.089024 22.57235 0.03428947 0.6043521 1.286215 0.002761350
## 115 8.113638 22.57432 0.03426871 0.6066691 1.283167 0.002631344
## 116 8.104102 22.57147 0.03421956 0.6047509 1.263558 0.002643695
## 117 8.085777 22.56867 0.03422197 0.6020866 1.259511 0.002719788
## 118 8.117772 22.57388 0.03423595 0.6066132 1.271117 0.002596183
## 119 7.632133 22.35861 0.03419148 0.5103708 1.262148 0.005699144
## 120 7.653227 22.40612 0.03418587 0.5289110 1.274611 0.005470742
## 121 7.679807 22.44377 0.03418710 0.5442404 1.285421 0.005268313
## 122 7.752634 22.50373 0.03420713 0.5702464 1.307663 0.004889237
## 123 7.705641 22.47525 0.03420189 0.5573022 1.300771 0.005231432
## 124 7.741408 22.48616 0.03418348 0.5630566 1.290899 0.004708992
## 125 7.807873 22.52078 0.03419986 0.5788447 1.303710 0.004354054
## 126 7.849828 22.52960 0.03419239 0.5835414 1.298763 0.003989176
## 127 7.868375 22.52751 0.03416603 0.5829817 1.283490 0.003720643
## 128 7.799326 22.50012 0.03415602 0.5701413 1.276727 0.004101735
## 129 7.724586 22.45353 0.03414840 0.5497171 1.266164 0.004582154
## 130 7.659449 22.37149 0.03415465 0.5168152 1.249295 0.005174357
## 131 7.660959 22.13551 0.03423078 0.4369002 1.186394 0.005825057
## 132 7.651473 22.25756 0.03414807 0.4765709 1.212502 0.005372958
## 133 7.678869 22.33712 0.03411141 0.5047827 1.225082 0.004922599
## 134 7.699846 22.39284 0.03411612 0.5256746 1.239725 0.004698500
## 135 7.754788 22.45066 0.03411603 0.5490425 1.250514 0.004261825
## 136 7.826030 22.49753 0.03412624 0.5693305 1.260894 0.003829193
## 137 7.955586 22.54331 0.03415405 0.5912152 1.274749 0.003229711
## 138 8.045038 22.56124 0.03415371 0.6004439 1.270085 0.002838828
## 139 7.972993 22.54389 0.03413385 0.5913783 1.263685 0.003116726
## 140 7.872204 22.50448 0.03410274 0.5725042 1.250056 0.003556161
## 141 7.786189 22.44577 0.03407990 0.5470814 1.235753 0.004030086
## 142 7.689798 22.26880 0.03408051 0.4811016 1.201009 0.004944878
## 143 7.686581 22.20129 0.03409430 0.4589797 1.186664 0.005147477
## 144 7.674517 22.14539 0.03416206 0.4412189 1.180503 0.005489079
## 145 7.736932 22.10127 0.03401697 0.4291697 1.157246 0.005045749
## 146 7.721798 22.22687 0.03402383 0.4674897 1.182527 0.004797429
## 147 7.748726 22.17319 0.03395086 0.4500800 1.165865 0.004747578
## 148 7.756942 22.28014 0.03396764 0.4843813 1.185002 0.004450208
## 149 7.765398 22.34451 0.03399965 0.5071950 1.199810 0.004278136
## 150 7.782598 22.38838 0.03401557 0.5235715 1.209524 0.004107591
## 151 7.837521 22.45553 0.03403934 0.5503014 1.223293 0.003736264
## 152 7.810138 22.42499 0.03402527 0.5378287 1.216421 0.003909328
## 153 7.884668 22.49241 0.03406263 0.5663666 1.232952 0.003485749
## 154 7.961869 22.53145 0.03409504 0.5845788 1.244613 0.003141870
## 155 8.030935 22.55408 0.03412522 0.5961641 1.253608 0.002877740
## 156 7.996781 22.53502 0.03409690 0.5827827 1.231228 0.003042449
## 157 7.973655 22.52436 0.03407228 0.5790960 1.230094 0.003120216
## 158 7.914841 22.49566 0.03403834 0.5659879 1.222594 0.003361816
## 159 7.956961 22.52211 0.03407051 0.5789522 1.232689 0.003169844
## 160 7.863123 22.44674 0.03399686 0.5448186 1.209987 0.003642236
## 161 7.805867 22.35727 0.03394283 0.5100171 1.191879 0.004044266
## 162 7.792273 22.28136 0.03388122 0.4823322 1.175614 0.004260508
## 163 7.773525 22.18100 0.03387100 0.4512562 1.161356 0.004576837
## 164 7.785727 22.14108 0.03379997 0.4383433 1.152331 0.004605158
## 165 7.806756 21.98780 0.03370065 0.3971126 1.130530 0.004906458
## 166 7.789603 22.05382 0.03377486 0.4142226 1.140319 0.004812733
## 167 7.792594 21.98708 0.03387741 0.3981260 1.133215 0.005019732
## 168 7.791164 21.95220 0.03404172 0.3900929 1.130345 0.005201062
## 169 7.772801 21.95562 0.03419456 0.3907170 1.133503 0.005409689
## 170 7.808551 21.89437 0.03425847 0.3766921 1.122964 0.005405158
## 171 7.866580 21.79210 0.03476391 0.3554618 1.109002 0.005439864
## 172 7.859142 21.86732 0.03457951 0.3650025 1.111913 0.005471933
## 173 7.843784 21.93159 0.03306298 0.3789830 1.118372 0.004969276
## 174 7.833579 21.99139 0.03320454 0.3916947 1.124214 0.004839394
## 175 7.817765 21.98868 0.03353210 0.3963581 1.128785 0.004844933
## 176 7.797035 22.16837 0.03375133 0.4448977 1.153587 0.004485273
## 177 7.804541 22.15549 0.03368836 0.4395171 1.149523 0.004485843
## 178 7.831620 22.28881 0.03373601 0.4705682 1.160749 0.004174662
## 179 7.829171 22.30228 0.03378759 0.4813284 1.167529 0.004101587
## 180 7.818590 22.30090 0.03382702 0.4852519 1.171972 0.004116655
## 181 7.829444 22.36688 0.03390737 0.5107089 1.186904 0.003932379
## 182 7.848979 22.40493 0.03393607 0.5249840 1.193963 0.003783662
## 183 7.859978 22.40320 0.03391001 0.5192265 1.187114 0.003778291
## 184 7.887471 22.45311 0.03397517 0.5427952 1.201899 0.003563062
## 185 7.846575 22.36815 0.03386451 0.5037638 1.178558 0.003912338
## 186 7.861556 22.40101 0.03389828 0.5100268 1.180612 0.003847025
## 187 7.899205 22.46266 0.03398840 0.5415167 1.199318 0.003549628
## 188 7.963114 22.52529 0.03411099 0.5710694 1.220760 0.003229571
## 189 7.989816 22.54255 0.03417254 0.5792129 1.230909 0.003135338
## 190 8.027567 22.55894 0.03423229 0.5910943 1.247302 0.002983625
## 191 8.001190 22.55737 0.03426545 0.5853027 1.245815 0.003137264
## 192 7.987880 22.55189 0.03424003 0.5806142 1.237733 0.003185428
## 193 7.961214 22.55358 0.03428668 0.5732841 1.236144 0.003365227
## 194 7.971280 22.55651 0.03429687 0.5773746 1.241666 0.003319763
## 195 7.983426 22.55938 0.03430645 0.5818937 1.248413 0.003266258
## 196 7.959657 22.56068 0.03433695 0.5754291 1.245682 0.003425298
## 197 7.939404 22.56601 0.03438051 0.5705355 1.246972 0.003593287
## 198 7.934369 22.56258 0.03435642 0.5658230 1.236969 0.003585531
## 199 7.908225 22.57678 0.03442494 0.5570800 1.235231 0.003815186
## 200 7.916297 22.57394 0.03442302 0.5628109 1.243421 0.003777831
## 201 7.900965 22.58384 0.03446208 0.5576141 1.242159 0.003921235
## 202 7.888318 22.58201 0.03441622 0.5417082 1.216503 0.003925962
## 203 7.937961 22.54563 0.03426097 0.5628582 1.224445 0.003481751
## 204 7.923497 22.54251 0.03425706 0.5555691 1.218250 0.003571040
## 205 7.904305 22.54390 0.03427183 0.5453128 1.211238 0.003713858
## 206 7.914347 22.52470 0.03418830 0.5487161 1.209428 0.003591499
## 207 7.933082 22.52906 0.03418508 0.5584723 1.215879 0.003461836
## 208 7.930186 22.51392 0.03411962 0.5563099 1.211089 0.003438151
## 209 7.911291 22.50716 0.03412266 0.5458899 1.204733 0.003573399
## 210 7.903700 22.48764 0.03406267 0.5410022 1.199791 0.003591100
## 211 7.916888 22.49026 0.03404672 0.5500432 1.204783 0.003473231
## 212 7.891354 22.46443 0.03400794 0.5332195 1.194083 0.003649609
## 213 7.882057 22.47205 0.03404642 0.5247191 1.190171 0.003772457
## 214 7.895865 22.50558 0.03413633 0.5361953 1.199326 0.003703018
## 215 7.880553 22.49766 0.03412427 0.5242472 1.191720 0.003835326
## 216 7.879475 22.51506 0.03417884 0.5248077 1.193489 0.003872016
## 217 7.868055 22.52000 0.03419372 0.5151887 1.187526 0.003995973
## 218 7.893325 22.52821 0.03421687 0.5363120 1.202300 0.003768700
## 219 7.889490 22.54843 0.03428925 0.5361260 1.205060 0.003838316
## 220 7.876956 22.56123 0.03432499 0.5277920 1.199850 0.003965051
## 221 7.881036 22.57042 0.03436033 0.5326869 1.205137 0.003947358
## 222 7.880083 22.58232 0.03440210 0.5341203 1.208231 0.003978777
## 223 7.878037 22.59544 0.03444757 0.5352177 1.211543 0.004024111
## 224 7.869626 22.60217 0.03444783 0.5272041 1.203343 0.004093340
## 225 7.862232 22.62595 0.03450031 0.5223417 1.199889 0.004189178
## 226 7.861231 22.61683 0.03446710 0.5192805 1.196015 0.004183863
## 227 7.868191 22.59333 0.03441446 0.5239683 1.199235 0.004089892
## 228 7.866725 22.58119 0.03437273 0.5205301 1.195139 0.004085801
## 229 7.859627 22.60855 0.03443608 0.5158647 1.192185 0.004188650
## 230 7.857624 22.59477 0.03438993 0.5114710 1.187656 0.004195564
## 231 7.863760 22.56645 0.03432317 0.5154725 1.190038 0.004097706
## 232 7.859761 22.54921 0.03426777 0.5091199 1.184362 0.004124324
## 233 7.853053 22.57486 0.03432266 0.5031444 1.180221 0.004235481
## 234 7.855348 22.52729 0.03420340 0.5017767 1.178499 0.004160674
## 235 7.849427 22.53101 0.03420287 0.4940456 1.172961 0.004253309
## 236 7.845176 22.53380 0.03419815 0.4870259 1.167791 0.004338469
## 237 7.845175 22.48629 0.03408897 0.4836174 1.165860 0.004299054
## 238 7.841309 22.48797 0.03408344 0.4759427 1.160519 0.004396808
## 239 7.838649 22.46827 0.03403191 0.4669653 1.154627 0.004482811
## 240 7.836543 22.37424 0.03383777 0.4561555 1.150001 0.004471674
## 241 7.839996 22.39206 0.03388178 0.4688298 1.157578 0.004327426
## 242 7.838619 22.39166 0.03387823 0.4648342 1.155126 0.004380240
## 243 7.842695 22.41469 0.03393138 0.4755919 1.161407 0.004277897
## 244 7.848716 22.44855 0.03400812 0.4873967 1.168415 0.004187633
## 245 7.859202 22.50793 0.03415734 0.5047449 1.180185 0.004090664
## 246 7.856202 22.46207 0.03404120 0.4985319 1.175209 0.004076288
## 247 7.863995 22.49094 0.03411330 0.5088577 1.182199 0.004007893
## 248 7.866452 22.45935 0.03402903 0.5099512 1.181627 0.003932018
## 249 7.868484 22.42736 0.03394688 0.5133738 1.182480 0.003836181
## 250 7.858320 22.41510 0.03393003 0.5010205 1.175700 0.003964571
## 251 7.848016 22.40159 0.03390623 0.4850557 1.166950 0.004137096
## 252 7.848557 22.37302 0.03385191 0.4892265 1.169231 0.004045120
## 253 7.840568 22.33606 0.03378725 0.4760341 1.162407 0.004151680
## 254 7.838923 22.33535 0.03377846 0.4689361 1.158444 0.004235028
## 255 7.835407 22.30628 0.03371785 0.4573347 1.152366 0.004340101
## 256 7.832624 22.27553 0.03368248 0.4583809 1.153767 0.004285703
## 257 7.829941 22.22905 0.03358752 0.4424943 1.145872 0.004423023
## 258 7.828547 22.11996 0.03338694 0.4174674 1.134812 0.004616247
## 259 7.832074 22.15842 0.03342443 0.4208610 1.135192 0.004627755
## 260 7.833672 22.27813 0.03365076 0.4418997 1.143878 0.004511631
## 261 7.835323 22.33171 0.03375422 0.4412458 1.141567 0.004625303
## 262 7.836413 22.14432 0.03341823 0.4114027 1.129561 0.004778491
## 263 7.846404 21.96175 0.03320305 0.3815350 1.118340 0.004999533
## 264 7.853743 21.89261 0.03347801 0.3703233 1.114149 0.005127411
## 265 7.849100 22.07150 0.03406740 0.3902783 1.118216 0.005219842
## 266 7.851934 21.96995 0.03390643 0.3788893 1.115846 0.005190184
## 267 7.847125 22.06015 0.03372091 0.3909130 1.119397 0.005096435
## 268 7.843016 22.06656 0.03344136 0.3950450 1.122110 0.004966332
## 269 7.841939 22.19001 0.03372092 0.4086669 1.125037 0.004987385
## 270 7.836613 22.44008 0.03396916 0.4490226 1.142308 0.004704216
## 271 7.836817 22.46904 0.03402425 0.4562632 1.146587 0.004639420
## 272 7.837992 22.58400 0.03425620 0.4626405 1.145928 0.004743841
## 273 7.839394 22.72597 0.03455146 0.4749124 1.148793 0.004790208
## 274 7.840241 22.75572 0.03460197 0.4741876 1.145772 0.004853913
## 275 7.838725 22.53730 0.03417777 0.4514890 1.138767 0.004855336
## 276 7.841427 22.75973 0.03459592 0.4704772 1.141254 0.004933471
## 277 7.839955 22.56403 0.03424377 0.4510015 1.136393 0.004926593
## 278 7.839077 22.41701 0.03399057 0.4366406 1.133502 0.004901061
## 279 7.840395 22.43139 0.03405732 0.4355044 1.131518 0.004967986
## 280 7.841807 22.60291 0.03434162 0.4514146 1.133779 0.005016430
## 281 7.843185 22.79559 0.03465209 0.4704526 1.138014 0.005016240
## 282 7.844780 22.84867 0.03473190 0.4734920 1.136573 0.005074184
## 283 7.848027 22.95524 0.03484206 0.4801620 1.133899 0.005201850
## 284 7.846005 22.86260 0.03474070 0.4729698 1.134482 0.005135862
## 285 7.844285 22.61106 0.03442899 0.4482076 1.129819 0.005147652
## 286 7.842923 22.41445 0.03415100 0.4298539 1.127686 0.005086427
## 287 7.841235 22.33706 0.03395345 0.4243514 1.128299 0.005000366
## 288 7.845247 22.26794 0.03412717 0.4124133 1.122958 0.005176622
## 289 7.845225 22.53971 0.03441577 0.4394239 1.126986 0.005221490
## 290 7.847382 22.19938 0.03425342 0.4034639 1.120346 0.005281017
## 291 7.845735 22.44007 0.03441132 0.4283241 1.124461 0.005294870
## 292 7.846625 22.73484 0.03462643 0.4575626 1.128949 0.005292429
## 293 7.849316 22.99612 0.03485511 0.4823921 1.132450 0.005299473
## 294 7.859441 23.53285 0.03541291 0.5310627 1.139382 0.005328000
## 295 7.853275 23.24564 0.03509975 0.5053556 1.135605 0.005311569
## 296 7.854093 23.30633 0.03526420 0.5117497 1.137955 0.005245243
## 297 7.862523 23.66248 0.03572908 0.5429834 1.142724 0.005280101
## 298 7.851447 23.20388 0.03521234 0.5037250 1.138621 0.005180496
## 299 7.854282 23.39391 0.03569463 0.5220514 1.144919 0.005135704
## 300 7.849487 23.18498 0.03535982 0.5047348 1.143373 0.005083215
## 301 7.844930 22.99012 0.03503580 0.4900521 1.143859 0.004996325
## 302 7.847920 23.17966 0.03545148 0.5070740 1.148070 0.005015191
## 303 7.853316 23.42336 0.03595043 0.5270624 1.151429 0.005073021
## 304 7.849197 23.28311 0.03573609 0.5174565 1.153316 0.005004013
## 305 7.843119 22.96701 0.03505471 0.4923857 1.149616 0.004906330
## 306 7.840869 22.87693 0.03490200 0.4907430 1.155599 0.004794061
## 307 7.843889 23.05764 0.03528865 0.5020751 1.154606 0.004897618
## 308 7.842626 23.01841 0.03522924 0.5021153 1.158738 0.004844903
## 309 7.844830 23.16593 0.03556033 0.5133926 1.161600 0.004904930
## 310 7.849145 23.38391 0.03598887 0.5293461 1.164884 0.005002989
## 311 7.857571 23.59865 0.03631215 0.5416790 1.155756 0.005107934
## 312 7.857126 23.60571 0.03635065 0.5430756 1.158944 0.005100759
## 313 7.865422 23.82107 0.03652719 0.5584983 1.155764 0.005190879
## 314 7.861639 23.69360 0.03629791 0.5478347 1.151483 0.005167605
## 315 7.864301 23.75011 0.03616374 0.5514303 1.148094 0.005222015
## 316 7.871658 23.95846 0.03624833 0.5683509 1.149428 0.005293699
## 317 7.872282 23.98675 0.03650080 0.5710335 1.153696 0.005266803
## 318 7.877267 24.10297 0.03645211 0.5802387 1.154898 0.005323694
## 319 7.873023 23.98403 0.03595756 0.5700353 1.147292 0.005366119
## 320 7.880031 24.16256 0.03602896 0.5849226 1.154707 0.005450761
## 321 7.880208 24.16567 0.03626352 0.5851909 1.155238 0.005389087
## 322 7.878515 24.13434 0.03643356 0.5828093 1.159543 0.005357183
## 323 7.875340 24.07010 0.03653669 0.5778900 1.160095 0.005311348
## 324 7.858812 23.68634 0.03643282 0.5498852 1.163496 0.005137548
## 325 7.853963 23.55927 0.03626093 0.5414562 1.165723 0.005085671
## 326 7.849877 23.45723 0.03607963 0.5356336 1.169539 0.005046746
## 327 7.853089 23.57796 0.03621396 0.5438708 1.170795 0.005118734
## 328 7.850066 23.49306 0.03609222 0.5386386 1.172222 0.005080644
## 329 7.850183 23.53056 0.03608211 0.5416311 1.175039 0.005126897
## 330 7.854569 23.63376 0.03623732 0.5477880 1.171940 0.005161334
## 331 7.858586 23.72567 0.03635051 0.5536884 1.169754 0.005193658
## 332 7.864048 23.85836 0.03637097 0.5629346 1.169570 0.005273677
## 333 7.872717 24.02647 0.03645536 0.5747938 1.164731 0.005326205
## 334 7.871773 24.01291 0.03637819 0.5738708 1.166472 0.005345860
## 335 7.874609 24.06693 0.03632746 0.5778195 1.165215 0.005380395
## 336 7.867984 23.95128 0.03623625 0.5694950 1.170142 0.005373556
## 337 7.865493 23.89548 0.03632924 0.5655566 1.170103 0.005311162
## 338 7.861617 23.83224 0.03619682 0.5613417 1.173491 0.005329320
## 339 7.858048 23.74872 0.03620983 0.5557476 1.174068 0.005267054
## 340 7.849563 23.54560 0.03603481 0.5430289 1.177696 0.005168716
## 341 7.852986 23.63525 0.03611775 0.5485106 1.176358 0.005216581
## 342 7.618733 22.31209 0.03431246 0.4828039 1.292801 0.007290257
## 343 7.620346 22.30741 0.03422251 0.4903649 1.256009 0.006116368
## 344 7.489675 22.42977 0.03438606 0.5110758 1.358705 0.008666885
## 345 7.446681 22.47494 0.03440033 0.5271089 1.366409 0.008842446
## 346 7.972693 22.56489 0.03435504 0.5928850 1.302107 0.003569201
## 347 7.634933 22.39057 0.03422082 0.5209876 1.288716 0.006116593
## 348 7.729173 22.35725 0.03406188 0.5126871 1.215550 0.004489801
## 349 7.815572 21.86922 0.03459762 0.3721966 1.122801 0.005639676
## 350 7.857996 21.80038 0.03459980 0.3578092 1.110787 0.005370024
## 351 7.859719 21.84592 0.03401293 0.3631377 1.111558 0.005234218
## 352 7.849913 22.67840 0.03467345 0.4493717 1.126798 0.005398281
## 353 7.849302 21.84226 0.03381187 0.3650018 1.113777 0.005184734
## 354 7.888251 22.59572 0.03448402 0.5483989 1.230508 0.004001169
## 355 7.880279 22.60183 0.03448332 0.5405022 1.219925 0.004040102
## 356 7.904643 22.57025 0.03438127 0.5507929 1.223185 0.003789128
## 357 7.872806 22.61788 0.03452575 0.5373113 1.218997 0.004130056
## 358 7.897121 22.56391 0.03435099 0.5441275 1.214734 0.003816325
## 359 7.874868 22.53920 0.03425291 0.5232162 1.194342 0.003952078
## 360 7.933279 22.56460 0.03433291 0.5622091 1.229492 0.003583631
## 361 7.853830 23.43986 0.03586860 0.5294313 1.152935 0.005082040
## 362 7.866847 23.89239 0.03652997 0.5648569 1.163614 0.005238604
## 363 7.585887 22.52606 0.03446174 0.5346649 1.345018 0.007793574
## 364 7.841383 22.33814 0.03380740 0.4873609 1.168981 0.004037820
## 365 7.468683 22.45532 0.03440032 0.5179380 1.362888 0.008753201
## 366 7.900712 22.56071 0.03437738 0.5861834 1.313313 0.004228143
## 367 7.868319 23.84531 0.03557715 0.5578915 1.144591 0.005440031
cor_5 <- rcorr(as.matrix(final))
M <- cor_5$r
p_mat <- cor_5$P
corrplot(M, type = "upper", order = "hclust", 
         p.mat = p_mat, sig.level = 0.01)
## Error in data.frame(..., check.names = FALSE): arguments imply differing number of rows: 78, 66

col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(M, method = "color", col = col(200),  
         type = "upper", insig = "blank", order = "hclust", # Add coefficient of correlation
         tl.col = "darkblue", tl.srt = 45, #Text label color and rotation
         # Combine with significance level
         p.mat = p_mat, sig.level = 0.01,  
         # hide correlation coefficient on the principal diagonal
         diag = FALSE 
         )

viol <- final %>%
 pivot_longer(
   cols = ends_with("per"),
   names_to = "ethnicity",
   values_to = "percentage",
   values_drop_na = TRUE
 )

viol <- subset(viol, select = -pop_count)
viol$ethnicity <- as.factor(viol$ethnicity)
library(ggplot2)
# Basic violin plot
p <- ggplot(viol, aes(x=ethnicity, y=percentage)) + 
  geom_violin(trim = TRUE)
p

corrplot(M, p.mat = p_mat, method = 'circle', type = 'lower', insig='blank',
         addCoef.col ='black', number.cex = 0.6, order = 'AOE', diag=FALSE)

# Get land area for tracts to calculate population density
sepabg <- tracts(state = "42", county = c("101", "017", "029", "045",
                                                "091"))
st_geometry(sepabg) <- NULL
head(sepabg)
##    STATEFP COUNTYFP TRACTCE       GEOID    NAME             NAMELSAD MTFCC
## 14      42      017  100302 42017100302 1003.02 Census Tract 1003.02 G5020
## 15      42      017  100304 42017100304 1003.04 Census Tract 1003.04 G5020
## 16      42      017  100306 42017100306 1003.06 Census Tract 1003.06 G5020
## 17      42      017  100307 42017100307 1003.07 Census Tract 1003.07 G5020
## 18      42      017  100401 42017100401 1004.01 Census Tract 1004.01 G5020
## 19      42      017  100402 42017100402 1004.02 Census Tract 1004.02 G5020
##    FUNCSTAT   ALAND  AWATER    INTPTLAT     INTPTLON
## 14        S 6861295 2033655 +40.0847565 -074.8915551
## 15        S 4018834   99694 +40.1122065 -074.8959059
## 16        S 1571083       0 +40.1061893 -074.8836008
## 17        S 2314482   71823 +40.0984599 -074.8984784
## 18        S 3245952    6352 +40.1348354 -074.8731389
## 19        S 3412871   20737 +40.1380006 -074.8558212
sepabg$A_KM2 <- sepabg$ALAND*1e-6

sf2_raw <- get_acs(geography = "tract", year = 2019, 
                         variables = c("B19001_001",  # Total households
                                       "B19001_002",   # Total <$10,000
                                       "B19001_006",   # Total 25000-29999
                                       "B19001_010",   # Total 45000-49999
                                       "B19001_013",   # Total 75000-99999
                                       "B19001_015",  # Total 125000-149999
                                       "B19001_017"   # Total >200000
                                       ),  
                         state = "PA", county = c("Philadelphia", "Bucks", 
                                                  "Chester", "Delaware", 
                                                  "Montgomery"))
## Getting data from the 2015-2019 5-year ACS
sf2 <- sf2_raw %>%
  group_by(GEOID) %>%
  summarise(house_tot = estimate[variable == "B19001_001"],
            pov_count = estimate[variable == "B19001_002"],
            low_count = estimate[variable == "B19001_006"],
            lower_med_count = estimate[variable == "B19001_010"],
            upper_med_count = estimate[variable == "B19001_013"],
            upper_count = estimate[variable == "B19001_015"],
            v_upper_count = estimate[variable == "B19001_017"]
            ) %>%
  mutate(pov_per = pov_count/house_tot*100,
         low_per = low_count/house_tot*100,
         lower_med_per = lower_med_count/house_tot*100,
         upper_med_per = upper_med_count/house_tot*100,
         upper_per = upper_count/house_tot*100,
         v_upper_per = v_upper_count/house_tot*100) %>%
  dplyr::select(GEOID, pov_per, low_per, lower_med_per, upper_med_per, upper_per, v_upper_per)
final <- merge(sf2, census_final, by = "GEOID")
final[] <- lapply(final, function(x) {
    if(is.factor(x)) as.numeric(as.character(x)) else x
})
sapply(final, class)
##         GEOID       pov_per       low_per lower_med_per upper_med_per 
##   "character"     "numeric"     "numeric"     "numeric"     "numeric" 
##     upper_per   v_upper_per          pm25          pm10           no2 
##     "numeric"     "numeric"     "numeric"     "numeric"     "numeric" 
##            o3            co           so2            pb 
##     "numeric"     "numeric"     "numeric"     "numeric"
final <- final %>% dplyr::select(-GEOID, -pm10)
x <- cor(final, use = "complete.obs")
col <- colorRampPalette(c("darkorange", "white", "steelblue"))(20)
corrplot(x, type = "upper", order = "hclust", col = col)